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ABSTRACT 

We present a detailed analysis of the local evolution of 206 Lagrangian Volumes (LVs) se¬ 
lected at high redshift around galaxy seeds, identified in a large-volume A cold dark matter 
(ACDM) hydrodynamical simulation. The LVs have a mass range of 1 — 1500 x 
We follow the dynamical evolution of the density field inside these initially spherical LVs 
from z = 10 up to ziow = 0.05, witnessing highly non-linear, anisotropic mass rearrange¬ 
ments within them, leading to the emergence of the local cosmic web (CW). These mass 
arrangements have been analysed in terms of the reduced inertia tensor I^j, focusing on the 
evolution of the principal axes of inertia and their corresponding eigendirections, and paying 
particular attention to the times when the evolution of these two structural elements declines. 
In addition, mass and component effects along this process have also been investigated. We 
have found that deformations are led by dark matter dynamics and they transform most of 
the initially spherical LVs into prolate shapes, i.e. hlamentary structures. An analysis of the 
individual freezing-out time distributions for shapes and eigendirections shows that first most 
of the LVs hx their three axes of symmetry (like a skeleton) early on, while accretion flows 
towards them still continue. Very remarkably, we have found that more massive LVs fix their 
skeleton earlier on than less massive ones. We briefly discuss the astrophysical implications 
our hndings could have, including the galaxy mass-morphology relation and the effects on the 
galaxy-galaxy merger parameter space, among others. 

Key words: gravitation - hydrodynamics - methods: numerical - galaxies: formation - cos¬ 
mology: theory - large-scale structure of Universe 


1 INTRODUCTION 


Over the last few decades, galaxy surveys such as the Two-degree- 
Field Galaxy Redshift Survey (2dFGRS; Colless^^e£^alJ|200n, the 
Sloan Digital Sky Survey (SDSS; e.g. Tegmark et al. 2004^, the 


Two-Micron All-Sky Survey (2MASS;|Huchra e t al.|2005| )md the 
6dFGS ( [Jones et al.|2004 ^ have revealed that galaxies gather in an 
intricate network, the so-called cosmic web (CW, after jBond, Kof-j 
[man, & Pogosyan||1996[ , made of filaments, walls, nodes which 
surround vast empty regions, the voids ( |Zel’dovi^ |1970[ jShan-j 
jdarin & Zeldovich|1989| l. These structures can be found on scales 
from a few to hundreds of megaparsecs and include huge flat struc¬ 
tures like the Great Wall ( [Geller & Huchra||1989| ( and the SDSS 
Great Wall ( jGott et al.|2005| (, the largest known structure in the lo¬ 
cal Universe, with a size larger than 400/i“^ Mpc and enormous 
empty regions like the Bootes void [Kirshner et al.|19^|1987[ . 

These results have been complemented by mappings of the 
dark matter (DM) spatial distribution through weak lensing obser- 
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vations like the Hubble Space Telescope Cosmic Evolution Sur¬ 
vey (COSMOS; [Massey et al.||2007l > and recent results from the 
Canada-France-Hawaii Telescope Lensing Survey (CFHTLenS; 
[Van Waerbeke et al.[2013[ l. 


Summing up, analyses of the current large scale distribution 
of galaxies and mass show that both are hierarchically organised 
into a highly interconnected network, displaying a wealth of struc¬ 
tures and substructures over a huge range of densities and scales. 
This web can be understood as the main feature of the anisotropi- 
cal nature of gravitational collapse [Peebles[1980[ , as well as of its 
intrinsic hierarchical character, and in fact it is the main dynamical 
engine responsible for structure formation in the Universe [Sheth [ 
[2004[[Sheth & van de Weygaert[20d4{ Shen et al.[2006l l, including 
galaxy scales (Dommguez-Tenreiro et al.[2011[(. 


According to the standard model of cosmology, large-scale 
structures observed in the Universe today are seeded by infinites¬ 
imal primordial density and velocity perturbations. The physical 
processes underlying their dynamical development until the CW 
emergence can be explained by theories and models on the gravita¬ 
tional instability, later on corroborated by a profusion of cosmolog- 
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ical simulations, the first of them purely A'^-body simulations (see 


e.g.,|Yepes, Dommguez-Tenreiro & Couchman|1992| Jenkins et al. 

19981 Pogosyan et al.|19981 Colberg, Krughoff & Connolly|2005[ 

Springel et al.|2005| Dolag et al.|2006), while recent ones include 

baryons and stellar physics too (see e.g.,|Dommguez-Tenreiro et al. 
20111|Metukietal.|2015J. 

Indeed, the advanced non-linear stages of gravitational insta- 
bility are described by the Adhesion Model (AM; see |Gurbatov 

& Saichev|1984[lGurbatov, Saichev & Shandarin|1989[|Shandarin 

& Zeldovich|19891|Gurbatov, Malakhov & Saichev|19921 Vergas- 

sola et al.ll 19941 and [Gurbatov, Saichev & Shandarin||20121 for a 

recent review), an extension of the popular non-linear Zeldovich 
Approximation (hereafter ZA; see|Zerdovich|1970|(. In comoving 
coordinates the ZA can be expressed as a mapping from the La¬ 
grangian space (the space of initial conditions g) into the Eulerian 
space (real space) described as a translation by a generalised irro- 
tational velocity-like vector (the displacement field s{q)) times the 
linear density growth factor D+{t), where the displacement can 
be written as a scalar potential gradient s{q) = —Vq4'(g). This 
approximation allows us to predict where singularities (locations 
with infinite density) will appear as cosmic evolution proceeds (i.e., 
the q points where the map has a vanishing determinant of the Ja¬ 
cobian matrix) and how they evolve into a sequence of caustics 
in real space. In this way, the ZA correctly but roughly describes 
the emergence of multistream flow regions, caustics and the struc- 

tural skeleton of the CW (jDoroshkevich, Ryaben’kii, & Shandarin 

T?73l IBuchertll 19891 fT99^handarin & ZeldovichlIT98'9l [Gales, 

Melott, & Shandarinll 19931 Melott, Pellman, & Shandarin|l 19941 

Melott, Shandarin, & Weinberg||1994| [Melott, Buchert, & Weib 

19951 ISahni & Coles|| 19951 [Yoshisato, Matsubara, & Morikawa 

19981|Yoshisato, Morikawa, Gouda, & Mouri|2006J. 

It is well known, however, that the ZA is not applicable once 
a substantial fraction of the mass elements are contained in multi¬ 
stream regions, because it predicts that caustics thicken and vanish 
due to multistreaming soon after their formation. One way of over¬ 
coming this issue is to introduce a small diffusion term in Zeldovich 
momentum equation, in such a way that it has an effect only when 
and where particle crossings are about to take place. This can be 
accomplished by introducing a non-zero viscosity, 0 , and then tak¬ 
ing the limit tz —>• 0: this is the AM, whose main advantage is that 
the momentum equation looks like the Burgers’ equation (|Burg- 
ers|l974| in the same limit, and hence its analytical solutions are 
known. A physically motivated derivation of the AM can be found 

in Buchert & Dominguez (|1998|(; Buchert, Dominguez & Perez- 


|Mercader| ( |1999| l; |Buchert & Dominguez| ( |2005| r 


The AM implies that, at a given scale, walls, filaments and 
nodes (i.e., the cosmic web elements) are successively formed, and 
then they vanish due to mass piling-up around nodes, to where mass 
elements travel through walls and filament^ Meanwhile, the same 
web elements emerge at larger and larger scales, and are erased at 
these scales after some time. Therefore, the AM conveniently de¬ 
scribes both the anisotropic nature of gravitational collapse and the 
hierarchical nature of the process. In addition, the AM indicates 
that the advanced stages of non-linear evolution act as a kind of 
smoothing procedure on different scales, by wiping mass accumu¬ 
lations off walls and filaments, first at small scales and later on at 
successively larger ones, to the advantage of nodes. Another impli- 


^ Recently confirmed in detail through CW element identification in large 
volume A-body simulations by|Cautun et al.||2014|. 


cation of the AM is that node centres (protohaloes at high z) lie on 
the former filaments at any 2 . 

A very interesting achievement of the AM is that the first suc¬ 
cessful reduction of the cosmic large scale structure to a geometri¬ 
cal skeleton was done in this approximation |Gurb atoy^Sajchey_& 


Shandarin|]^9^ Kofman, Pogosyan & Shandarin|]^9W Gurbatov, 
Saichev & ShandarinD2012|), see al^Hidding, Shandarin & 


de Weygaert pOM 1 . Later on Novikov, Colombi & Dore|P006[); 


Sousbie et aL| ( |2008a|b^ ; Sousbie, Colombi & Pichon|P009t ^ 


bI^P011|l;|Sousbie, Pichon & Kawahara| pOll^; |Arag6n-Calvo, 


van de Weygaert & Jones|P010[l and |Arag6n-Calvo et al.H2010[l 


also discussed the skeleton or spine of large-scale shuctures from 
purely topological constructions in a given density field. 

Recently, a growing interest to identify and analyse elements 
of the CW in A-body simulations, as well as in galaxy catalogues, 


et al.l2005 lAragon-Calvo et al.|2007alb 

2010 Hahnetal.l2007b|al 

Platen, van de Weygaert & Jones||2007 

Stoica, Martinez & Saar 

20071 lEorero-Romero et al.||20()‘)l |Wu, Batuski &. KhalilllSOOOl 

Aragon-Calvo, van de Weygaert & Jones||20101 IBond, Strauss & 

Cen|12010albl [Genovese et al. 

2 OI 2 I 

Tonzalez & Padilla[[20101 

Jones, van de Weygaert & Arag6n-Calvo[[2010[ [Stoica, Martinez 

& Saar[[2010[ [Hoffman et al.[ 

2012[ [Cautun, van de Weygaert & 

Jonesl[2013[ Tempel et al.[[20 

These methods and algorithms 


are motivated by the study of the influence of large scale struc¬ 
tures on galaxy formation (|Altay, Colberg & Croft|2006 |Arag6n- 


Calvortay200^ Halme£ay2007bJ^ Paz, Stasyszyn & Padilla 


2008r Hahn et al. [2009[ [Zhang et al. [2009 

[Codlowski, Panko & 

Flinl201 l|[Codis et al.|2012llLibeskind et a 

l.|2012 2013[ Aragon- 


Calvo^& ifang j2Ql^ Metuki et al.|2015| . In a recent paper, Cau 


tun et al. ( 2014^ have investigated the evolution of the CW from 


cosmological simulations, focusing on the global evolution of their 
morphological components and their halo content. 

From a dynamical point of view, |Hidding, Shandarin & van] 
|de Weygaert|p014] > go a step further by establishing the link be¬ 
tween the skeleton or spine of the CW, as described by the previous 
methods, and the development of the density field. In fact, they de¬ 
scribe for the first time the details of caustic emergence as cosmic 
evolution proceeds. Their main result is to show that all dynamical 
processes related to caustics happen at locations placed near a set of 
critical lines in Lagrangian space, that, when projected onto the Eu- 
lerian space, imply an increasing degree of connectedness among 
initially disjoint mass accumulations in walls or filaments, until a 
percolated structure forms, i.e., the spine or skeleton of the large 
scale mass distribution. These authors compare their results with 
two dimensional A-body simulations. Note that, due to the com¬ 
plexity of the problem, they first work in two dimensional spaces, 
where caustic emergence and percolation are described. Neverthe¬ 
less, they expect no important qualitative differences when three- 
dimensional spaces are considered instead. 

As we can see, in the last years different methods to quantify 
the cosmic web structure, classify its elements and study its emer¬ 
gence and evolution have been developed and applied. However, 
a detailed analysis of the local development of the density field 
around galaxy hosting haloes is still missing. This is of major im¬ 
portance because of its close connection to the problem of galaxy 
formation, in which case the effects of including gas processes need 
to be considered too. It is worth noting that neither the ZA nor the 
AM include gas effects in their description of CW dynamics. 

This analysis should first answer to the simplest questions re¬ 
lated to local shape deformation and spine emergence and the ori¬ 
entation of its main directions or symmetry axes around galaxy-to- 
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be objects. Besides, the very nature of these local processes, there 
are other interesting, simple, not-yet-elucidated related issues. For 
instance the characterisation of the times when deformation stops 
and orientation gets frozen, whether or not this local weh evolution 
is mass dependent (i.e., the mass of the halo-to-he) or not, and if 
different components (DM, hot gas, cold baryons) evolve in a sim¬ 
ilar way or there is a component segregation. We do not have at 
our disposal an analytical tool to perform such analyses, in conse¬ 
quence we need to resort to numerical simulations. 

In order to answer these questions, in this paper we investigate 
the impact of the local features of the Hubble flow imprinted on the 
deformation of initially spherical Lagrangian volumes (LVs) and 
the spine emergence, from high to low redshift. As known from 
previous studies, the local Huhble flow is neither homogeneous 
nor isotropic, on the contrary, it contains shear terms (and small- 
scale vorticity at its most advanced stages) that distort cosmolog¬ 
ical structures. We use cosmological hydrodynamical simulations 
to study the deformations of a sample of LVs through their reduced 
inertia tensor at different redshifts, which allows us to describe in 
a quantitative way the LV shape deformation and evolution, along 
with that of their symmetry axes. We analyse every component sep¬ 
arately, that is, we compute the reduced inertia tensor for DM, cold 
and hot haryons. 

This paper is organised as follows. In ^ we outline the simu¬ 
lation method and the algorithms used to study the deformations of 
LVs. A brief summary on the ZA, the CW emergence in 2D and the 
AM is given in ^ where some of their implications, useful in this 
paper, are also addressed. Some relevant details of the highly non¬ 
linear stages of gravitational instability, beyond the ZA or the AM 
are summarised in to help to understand how our results about 
the LV evolution can be explained in the light of these models. In 
the LV evolution is investigated in terms of the reduced inertia 
tensor eigenvectors, delaying the analysis in terms of its eigenval¬ 
ues to the next section, focused on the mass and component ef¬ 
fects and on the shape evolution of the selected LVs. In Sj^we study 
the freezing-out of eigendirections and shapes, presenting the dis¬ 
tribution of the corresponding freezing-out times and looking for 
mass effects. Possible scale effects on the previous results are dis¬ 
cussed in !|^ Finally, we present our summary, conclusions and 
discussion in 


2 SIMULATIONS AND METHODS 
2.1 Simulations 


The simulations analysed here have been run under the GALFOBS 
I and II projects. The GALFOBS (Galaxy Formation at Different 
Epochs and in Different Environments: Comparison with Observa¬ 
tional Data) project aims to study the generic statistical properties 
of galaxies in various environments and at different cosmological 
epochs. This project was a DEISA Extreme Computing Initiative 
(DEClfJ GALEOBS I was run at LRZ (Leibniz-Rechenzentrum) 
Munich, as a European project. Its continuation, GALEOBS 11, was 
run at the Barcelona Supercomputing Centre, Spain. 

All the runs were performed using P-DEVA, the parallelised 
of the DEVA code ([Serna, Domlnguez-Tenreiro & Saiz] 


20031. DEVA is an hybrid AP^M Lagrangian code, implemented 


with a multistep algorithm and smoothed particle hydrodynamics 
(SPH). The SPH version included in P-DEVA ensures energy and 
entropy conservation and, at the same time, guarantees a good de¬ 
scription of the forces and angular momentum conservation. How¬ 
ever, this advantage implies a gain in accuracy and an additional 
computational cost. Star formation (SE) is implemented through a 
Kennicutt-Schmidt-like law with a given density threshold, p,, and 
star formation efficiency c, ( [Martinez- Serrano et al.|200^ . 

The simulations have been carried out in the same periodic 
box of 80 Mpc side length, using 512® haryonic and 512® DM 
particles. Due to computational cost, these simulations only in¬ 
clude hydrodynamical calculation in a sub-box of 40 Mpc side. 
The evolution of matter follows the A cold dark matter (ACDM) 
model, with parameters Dm = 0.295, Db = 0.0476, Da = 0.705, 
h = 0.694, an initial power-law index n = 1, and erg = 0.852, 
taken from cosmic microwave background anisotropy dat^ (Dunk- 


[ley et al.|2009^ . The star formation parameters used were a density 
threshold pthres = 4.79 x 10“®®g cm~® and a star formation ef¬ 
ficiency c = 0.3. The mass resolution is rribar = 2.42 x Mq 
and rtiDM = 1.26 x 10®Mq and a spatial resolution of 1.1 kpc in 
hydrodynamical forces. More detailed information of these simu¬ 
lations can be found in |Onorbe et al.| ( |2011^ . 

It is noteworthy that no explicit feedback has been imple¬ 
mented in these simulations, but SE regulation through the values 
of the SF parameters. Nevertheless, the issues that will be discussed 
in this paper involve considerably larger characteristic scales than 
the ones related to stellar feedback. Therefore, it is unlike that the 
details of the star formation rate, and those of stellar feedback in 
particular, could substantially alter the conclusions of this paper. 


2.2 Methods 


We first describe how the LV sample around simulated galaxies has 
been built up. The first step is halo selection at ziqw = 0-05 by us- 


the SKID algorithrrPljWeinberg, Hernquist & Katz 


1997 I. This 


mg 

multi-step algorithm determines first the smoothed density field, 
then it moves particles upward along the gradient of this density 
field using a heuristic equation of motion that forces them to col¬ 
lect at local density maxima. Afterwards, it defines the approximate 
group to be the set of particles identified with an FOF algorithm 
with a linking length, b. Finally, particles not gravitationally bound 
to the groups identified in the previous step are removed. 

Specifically, we have selected a sample of 206 galaxy haloes 
from two runs of the GALFOBS simulations at Ziow, not involved 
in violent events at the halo scale at zio-w - Their virial radii rvir.iow 
and masses Mvir.iow at this redshift go from dwarf galaxies to 
galaxy groups, see the corresponding histograms in Fig.[^first row. 
The virial radius (rvir) is defined as the radius of the sphere enclos¬ 
ing an overdensity given by |Bryan & Norm^ ( |1998| l. 

Next, for each halo at ziow we have traced back all the particles 
inside the sphere defined by its respective rvir.iow to Zhigh = 10. 
Using the position of these particles at Zhigh we have calculated a 
new centre ry. Then, we have selected at Zbigh all the particles en¬ 
closed by a sphere of radius i?high = K x rvir,low, with K — 10 
around their respective centres fc (see first row of Fig. |^, and 
we have identified each of the DM and haryonic particles within 
these spherical volumes. These particles sample the mass elements 


^ The DEISA Extreme Computing Initiative was launched in May 2005 
by the DEISA Consortium, as a way to enhance its impact on science and 
technology 


® http://lambda.gsfc.nasa.gov/product/map/dr3/params/ 
lcdm_sz_lens_run_wmap5_bao_snall_lyapost.cfm 
^ http://www-hpcc.astro.washington.edu/tools/skid.html 
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Figure 1. Upper panels show the radius and mass distribution of the galaxy 
haloes at 2 io„ in our sample. Lower panels depict the same information for 
the selected LVs. 


whose deformations, stretchings, foldings, collapse and stickings 
we are to trace along cosmic evolution. They follow geodesic tra¬ 
jectories until they possibly get stuck and begin the formation of, or 
are accreted onto, a CW structure element. For this reason, we have 
termed them Lagrangian Volumes (LVs). It is worth noting at this 
point that we are following the evolution of individual LVs, each 
of them made of a fixed number of particles as they evolve. We do 
not trace the possible incorporation of off-LV mass elements that 
could happen along evolution as a consequence of mergers, infalls 
or other processes. Note also that, due to the very complex evolu¬ 
tion of the LVs, their borders are not well defined at 2 < 2high- 
Finally, a technical point to take into account is that the LVs should 
lie inside the hydrodynamical zoomed box. 

The choice K — 10 is motivated as a compromise between 
low K values, ensuring a higher number of LVs in the sample, and a 
high K, ensuring that LVs are large enough to meaningfully sample 
the CW emergence around forming galaxies. The possible effects 
that different K values could have in our results will be discussed 
in Section where we conclude that fT = 10 is the best choice 
among the three possibilities analysed. 

Afterwards, we have followed the dynamical evolution of 
these particles across different redshifts until they reach 2iow, i-e., 
we have followed the evolution (stretchings, deformations, fold¬ 
ings, collapse, stickings) of a set of 206 LVs from 2high until 2iow 
By construction, the mass of each of these sets of particles is con¬ 
stant across evolution, and its distribution is given in Fig.[T] second 
row, where we also show the distribution of their initial sizes at 

•2^high • 

The choice of initially spherically distributed sets of particles 
aims to unveil the anisotropic nature of the local cosmological evo¬ 
lution, illustrated in Fig.|^ where two examples of LVs at 2 = 10 
and their corresponding final shapes and orientations at 2iow are 
displayed. The mass of these LVs are 8.7 x (left-hand 

panels) and 4.4 x (right-hand panels), respectively. 

In this figure we note that, in both cases, a massive galaxy 


appears at 2iow in the central region of the LV. It turns out that, 
by construction, these galaxies are just those identified in the first 
step of the LV sample building-up, see above. We also notice that 
the LVs have evolved into a highly irregular mass organisation, in¬ 
cluding very dense subregions as well as other much less dense 
and even rarefied ones. Also, some changes of orientation of the 
emerging spines are visible, mainly in the lighter LV. In addition, 
the initial cold gaseous configuration at 2 = 10 has been trans¬ 
formed into a system where stars (in blue) appear at the densest 
subregions of the LVs. Hot gas (in red) particles are also present 
and constitute an important fraction of the LV mass (see sQfor an 
explanation about its origin). We also observe that the overall LV 
shape on the right-hand side of Fig. is highly elongated at 2iow 
and has a prolate-like or filamentary appearance, visually spanning 
a linear scale of ~ 9 Mpc long by 2 Mpc wide, while that on the 
left-hand side of Fig.|^still keeps a more wall-like structure. These 
shape transformations illustrate the highly anisotropic character of 
evolution under gravity. In this respect, it is worth mentioning that 
anisotropy is a generic property of gravitational collapse for non¬ 
isolated systems, as it was pointed out in early works by|Lin, Mestel| 
|& Shu|fT%5) ; [Ic^ ( [TW^ and |White & ^ ( fTW^ . 

As we mentioned in iJT] the deformation, stretching, folding, 
multistreaming and collapse of mass elements by cosmological 
evolution is predicted and described by the ZA, while AM adds 
a viscosity term making multistreaming regions to get stuck into 
dense configurations. In the following, we will introduce the math¬ 
ematical methods we use to quantify the local LV transformations 
illustrated in Fig.|^ 

To this end, we have calculated, at different redshifts, the re¬ 
duced inertia tensor of each LV relative to its centre of mass 






n = 1, ..., W 


( 1 ) 


where r„ is the distance of the n-th LV particle to the LV centre of 
mass and N is the total number of such particles. We have used this 
tensor instead of the usual one ( [Porciani, Dekel & Hoffman|2002a| l 
to minimise the effect of substructure in the outer part of the LV 
( |Gerhard|1983|[Bailin & Steinmetz|2005| l. In addition, the reduced 
inertia tensor is invariant under LV mass rearrangements in radial 
directions relative to the LV centre of mass. This property makes 
the Ilj tensor particularly suited to describe anisotropic mass de¬ 
formations as those predicted by the ZA and the AM and observed 
in Fig. 

In order to measure the LV shape evolution, first, we have cal¬ 
culated the principal axes of the inertia ellipsoid, a, b, and c, derived 
from the eigenvalues (Ai, with Ai A2 ^ A3) of the /jj tensor, so 
that ab'^ c (see|Gonzalez-GarcIa & van AlbadaH2005^), 


a = 



— Ai -I- As) 

2M 


b = 


5(A3 — A2 -f Ai) 
2M ’ 

5(Ai — As -f A2) 
2M ’ 


( 2 ) 


where M is the total mass of a given L\Q We denote the directions 
of the principal axes of inertia by e^, i = 1 , 2 ,3, where ei corre¬ 
spond to the major axis, 62 to the intermediate one and is to the 
minor axis. 

Afterwards, to quantify the deformation of these LVs, we have 


® Note that Ai + A 2 + A 3 = 2M and this implies = 5. 
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Figure 2. Left, shape evolution of a wall-like LV from 2 = 10 to 2iow = 0.05. Different columns are three projections of the same LV, with fixed axes 
taken oriented along the direction of the principal axes at 2iow Magenta points represent DM, green cold gas, red hot gas (T ^ 3 x 10“^ K) and blue stars. 
First row shows the initially spherical LV at 2 = 10, where DM and cold gas are represented in the same plot. Second, third, fourth and fifth group of panels 
illustrate the LV shape deformation across redshifts 2 = 3,1,0.5 and 0.05, where DM and baryonic components are split in different rows. Right, the same 
for a filament-like LV. The mass of the LVs are 8.7 x 10^^ Mq and 4.4 x 10^^ Mq, respectively. 
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computed the triaxiality parameter, T, ( |Franx, Illin g worth & de| 
|Zeeuw|1991|), defined as 


T = 


(1 - 

{l-c^/a?y 


(3) 


where T = Q corresponds to an oblate spheroid and T = 1 to 
a prolate one. An object with axis ratio cja > 0.9 has a nearly 
spheroidal shape, while one with c/a < 0.9 and T <0.3 has an 
oblate triaxial shape. On the other hand, an object with c/a < 0.9 
and T > 0.7 has a prolate triaxial shape (|Gonzalez-Garcfa et al.| 
[ 200 ^ . 

We have also calculated other parameters that measure shape 
deformation such as, ellipticity, e 


a’^ + + c^’ 


(4) 


that quantifies the deviation from sphericity, and prolateness, p 
a'^ + c^ - 2F 


that compares the prolateness versus the oblateness ^Bardeen_e£^ 


198^ Porciani, Dekel & Hoffman)2002b{|Springel, White & Hem 

quist|2004^ . In this case, a sphere has e = p = 0, a circular disc 
has e = 0.5, p = —0.5 and a thin filament has e = p = 1. Nearly 
spherical objects have e < 0.2 and |p| < 0.2. 

To sum up, we have performed the computation of the reduced 
inertia tensor, the principal axes of inertia, the eigendirections and 
the parameters T, e and p involving each of the selected LVs. Fur¬ 
thermore, we have repeated the same calculation for each compo¬ 
nent separately, viz. DM, cold and hot baryons. We consider hot 
gas as the particles shock heated to 3 x 10^ K. 


3 EVOLUTION UNDER THE ZA OR THE AM 


The advanced non-linear stages of gravitational instability are de- 


scribed by the adhesion model ^Gurbatov &. Saichev||1984 

Gur- 

jbatov, Saichev & Shandarin||I989||Shandarin & Zeldovich| 

00 


[Gurbatov, Malakhov & Saichev|1992||Vergassola et al.|1994^ , an 

extension of Zeldovich’s (1970) popular non-linear approximation. 
In this Section, we briefly revisit them as well as some of their im¬ 
plications, useful to understand the results that will be analysed in 
the next sections. 


3.1 The Zeldovich Approximation 

In comoving coordinates, Zeldovich’s approximation is given by 
the so-called Lagrangian map: 

Xi{q,t) = qi +D+{t)si{q), (6) 

where qi and Xi,i = 1,2,3 are comoving Lagrangian and 
Eulerian coordinates of fluid elements or particles sampling them, 
respectively (i.e., initial positions at time tin and positions at later 
times f); D+{t) is the linear density growth factor. As already men¬ 
tioned, it turns out that Si{q) can be expressed as the gradient of the 
displacement potential 'i’iq). 

The behaviour of D+{t) depends on the cosmological epoch. 
For the flat concordance cosmological model (see § 113’ at high 
enough z, when the Universe evolution is suitably described by 

the Einstein-de Sitter model, D+{t) = (3/5)(f/fi)^^®. Later on, 
,2 

when ~ 0 and the effects of the cosmological constant emerge 
(2 a — 0.684 orfA/fu = 0.554 for the cosmological model used in 


the simulations analysed here), D+ (t) is an exponential function of 
time. Finally, when the cosmological constant dominates, we have: 


D+{a{t)) oc ®45/6,2/3) 


Go \ 

Ga j 



o 1 G2 
\Im 

u^Ga 


(7) 


where 93a; is the incomplete /3 function. Go = 1 — Ga, Gm is the 
non-relativistic contribution to Go, and 


o'^Ga 


Go -I- o^Ga ’ 


(8) 


describing a frozen perturbation in the limit f —^ oo. 

Due to mass conservation, equation implies for the local 
density evolution: 


p(f, t) = 


Pb{t) 


[1 - D+{t)a{q)\[l - D+{t)P{^\[l - D+{t)'y{q)\ 


,(9) 


where r = a{t)x is the physical coordinate, pbit) the background 
density, and 'y{q) < f5{q) < a{q) are the eigenvalues of the lo¬ 
cal deformation tensor, dij{q) = — ■ Equationj^describes 

caustic formation in the ZA. Indeed, a caustic first appears when 
and where Dj,.{t)a(q) = 1 (i.e., a wall-like one), see details in § 
|3.2[ Mathematically, caustics at time t can be considered as singu¬ 
larities in the Lagrangian map (see equation]^ and more details in 
the next subsection). 


3.2 The CW Emergence in 2D 


The emergence of the cosmic skeleton as cosmic evolution pro¬ 
ceeds in the frame of the ZA is presented by [Kidding, Shandai mj 
|& van de Weygaert] p014| >. Due to the high complexity of the 
formalism involved, the authors restrict themselves to the two- 
dimensional equivalent of the ZA, providing us with the concepts, 
principles, language and processes needed as a first step towards a 
complete dynamical analysis of the CW emergence in the full three- 
dimensional space. In this subsection we give a brief summary of 
some of their results, useful to interpret some of our findings. 

In 2D, the complexity of the cosmic structure can be under¬ 
stood to a large extent from the properties of the a{q) landscape 
field, where a{q) is the largest eigenvalue of the deformation ten¬ 
sor dij{q), i,j = 1, 2. The role of the second eigenvalue j3{q) is 
much less relevant, except around the places where the haloes are 
to form. 

Of particular relevance are the A 3 lines in Lagrangian space, 
because they are the progenitors of the cosmic skeleton in Eule¬ 
rian space. Geometrically they can be defined as the locus of the 
points where the gradient of a (or /3) eigenvalue is normal to its 
corresponding eigenvector Ca (or ep). Alternatively, they can also 
be defined as the locus of the points where Ca (or ejg) is tangential 
to the contour level of the a{q) (or l3{q)) landscape field. 

The locations where collapse first occurs are around the max¬ 
ima of the a(q) field in Lagrangian space. These are the so-called 


singularities, after Arnold’s singularity classification (Arnold 


|1983| l. They are placed on the A 3 lines. Subsequently, the evolution 
under the ZA drives a gradual progression of Lagrangian collaps¬ 
ing regions, consisting, at a given time t, of those points such that 
a{q) = l/Z)+(f) or/3(^ = l/G+(f), according to the 2D version 
of equation]^ These isocontours lines are the so-called A^d) and 
A 2 (t) lines, and within them matter is multistreaming in Eulerian 
space, i.e., matter forms a fold caustic or pancake. 

The height of the a{q) landscape field portrays the collapse 
time for a local mass element. Indeed, at a given time t, points 
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where the A^f) and the A 3 lines meet, correspond to points in 
Eulerian space where a cusp singularity can be found (i.e., the tip of 
a caustic). The A^f) lines descend on the a{q) landscape field as 
time elapses, and in this way more and more mass elements get in¬ 
volved in the pancake. The pancake grows in Eulerian space, where 
the two cusp singularities at their tips move away from each other. 
A similar description can be made for the p{q) eigenvalue. 

Note that the height of either the A 2 (t) or the Aj (t) lines 
depends only on the D+{t) function, and not on the eigenvalue 
landscape fields. Therefore, the higher the a{q) landscape field, the 
earlier the corresponding pancake in the Eulerian space is formed. 
The same argument holds for the eigenvalue. 

Along the A 3 lines there are another types of extrema. Eirst, 
we have the Aj singularities or saddle points, after Arnold’s clas¬ 
sification. They are in-between two A;^ singularities and are local 
minima along the A 3 lines. They depict the places where two pan¬ 
cakes emerging from each of the A;^ points get connected, when 
the corresponding A 2 lines met the Aj singularities at their de¬ 
scent. This represents a first percolation event, and a first step to¬ 
wards the emergence of the CW spine. Eor the aforementioned rea¬ 
sons, the higher the a{q) landscape field, the earlier the percolation 
events will occur. 

The second type are the local maxima points q 4 , where the 
corresponding eigenvector is tangent to the A 3 lines, i.e., the so- 
called A4 singularities, or swallow tail according to |Amold| ( |1983| l. 
An A 4 singularity at q 4 exists only at a unique instant t 4 , when 
a(^) = l/D+(f 4 ). At this moment, the A 2 (f 4 ) line passes 
through A 4 , transforming the cusp singularity at the end of the Eu¬ 
lerian pancake into a swallow tail singularity. After that, there are 
three intersections of the A 2 (t) line with two A 3 lines, giving three 
connected cusp singularities in Eulerian space. Therefore, the A 4 
singularities are the connection points where disjoint pieces of A 3 
lines get connected in Eulerian space. Then, we get another perco¬ 
lation process. Once again, as explained above, the higher the a{q) 
landscape field, the earlier the percolation events will take place. 

This short summary illustrates some aspects of the effect that 
the height of the a{q) landscape field has on the time when simple 
percolation events occur in 2D, or, in a more general scope, when 
the CW spine emerges. The conclusion is simple: the higher the 
eigenvalue landscape, the earlier the percolation events take place. 
A similar effect can be expected in 3D, provided that the descrip¬ 
tion of the events connecting disjoint caustics in Eulerian space is 
not dramatically changed with respect to that in 2D. 

Pancake formation in Eulerian space entails an anisotropic 
mass rearrangement as matter flows normally to the a (or /3) pan¬ 
cake. These flows consist of mass elements within the A^f) (or 
Ajjf)) lines in Lagrangian space, and therefore they ideally do 
not stop while the A 2 lines keep on descending on the landscape. 
Similar ideas apply to other kind of caustic formation, implying 
shape transformations after the skeleton emergence. Note that mat¬ 
ter flows are predominantly anisotropic, except for the places where 
the haloes are to form, i.e. where flows become more isotropic. 


3.3 The Adhesion Model 

As it is well known, Zeldovich’s approximation is not applicable 
beyond particle crossing, because it predicts that caustics thicken 
and vanish due to multistreaming soon after their formation. How¬ 
ever, A-body simulations of large-scale structure formation indi¬ 
cate that long-lasting pancakes are indeed formed, near which par¬ 
ticles stick, i.e multistreaming did not take place. The adhesion 
model was formulated to incorporate this feature to Zeldovich’s ap¬ 


proximation, by introducing a small diffusion term in Zeldovich’s 
momentum equation, in such a way that it has an effect only when 
and where particle crossings are about to take place. This can be 
accomplished by introducing a non-zero viscosity, v, and then tak¬ 
ing the limit i/ ^ 0. This is the phenomenological derivation of 
the adhesion model. Physically motivated derivations can be found 
Buchert & Dommguez| jl998^, [Buchert, Dominguez & Perez 


in 


Mercader 


DommguezH20()5| >. 


1 1999|and others included in the review by Buchert & 


As in the Zeldovich approximation, in the adhesion model, 
the initial velocity field can be expressed as the gradient of a scalar 
potential field, 4>o((?), describing the spatial structure of the initial 
perturbation. It can be shown that the solutions for the velocity field 
behave just as those of Burgers’ equation ( |Burgers|1948|[T974) in 
the limit i/ —> 0 , whose analytical solutions are known. 

The most significant characteristic of Burgers’ equation so¬ 
lutions is that they are discontinuous and hence they unavoidably 
develop singularities, i.e., locations where at a given time the ve¬ 
locity field becomes discontinuous and certain particles coalesce 
into long-lasting very dense configurations with different geome¬ 
tries, i.e., caustics as in the ZA. The ideas explained in § |3.2| also 
apply here, but the main difference is that matter gets stuck forming 
very dense subvolumes (singularities) in Eulerian space, instead of 
forming multistreaming regions. In this way, a singularity occurs 
at the time t when a non-zero d-dimensional elemental volume V 
around a point q in the initial configuration is mapped to a d'- 
dimensional elemental volume around a point x{q, t) in Eulerian 
space with d' < d. In a three-dimensional space, these singulari¬ 
ties can be walls (with dimension d' = 2 ), filaments (d' = 1 ) and 
nodes (d' = 0 ). 

The AM model implies that, locally, walls are the first sin¬ 
gularities that appear, as denser small surfaces (the so-called pan¬ 
cakes). Later on, filaments form and grow until singularity per¬ 
colation and spine emergence ((Gurbatov, Saichev & Shandarin 
198^[Kofman, Pogosyan & Shandarm[1990||G^atov, Saichev & 
Shandarin]2012| . The singularity pattern implies the emergence of 
anisotropic mass flows towards the new formed singularities. Lo¬ 
cally, emerging walls are the first that attract flows from voids, then 
they host flows towards filaments, and, finally, filaments are the 
paths of mass towards nodes. In this way, at a given scale, walls 
and filaments tend to vanish as the mass piles up at nodes. In addi¬ 
tion, cells associated with the deepest minima of — tI>o(^. swallow 
up some of their neighbouring cells related to less deep minima, in¬ 
volving their constituent elements (i.e., walls, filaments and nodes), 
and causing their merging, as in the ZA. This is observed in sim¬ 
ulations as contractive flow deformations that erase subsbucture at 
small scales, as mentioned above, while the CW is still forming at 
larger scales. 

It is worth noting that Burgers’ equation solutions ensure the 
existence of regular points or mass elements at any time t, as those 
that have not yet been trapped into a caustic at t. Because of that, 
these regular mass elements are among the least dense in the den¬ 
sity distribution. Note, however, that due to the complex structure 
of the flow, singular (i.e., already trapped into a caustic) and regu¬ 
lar (i.e., not yet trapped) mass elements need not be spatially seg¬ 
regated, and in fact, they are mixed ideally at any scale. 


3.4 Further implications 

According to the ZA, we have 
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5Sp 


7q- • s = a(,j) + /3((j) + 7® = (tin). 


( 10 ) 


As suggested by the 2D analysis made in § |3.2[ the height of 
the a{q) landscape field in 3D portrays the collapse time for local 
mass elements (with a{q) the larger dij eigenvalue at q), as well as 
the time when different percolation events mark the emergence of 
the CW spine. Equation [T^ indicates that the eigenvalue landscape 
fields are closely related to the fluctuation field (FF) ^ at tin. 

It is well known that the number density of the FF peaks above 
a given threshold is considerably enhanced by the presence of a 
(positive) background field ( [Bardeen et al.|198^ , or, equivalently, 
when a large-scale varying field is added to —. Equation 


10 


tells 

us that such background would increase the height of the landscape 
fields, thereby speeding up percolation events responsible for the 
CW emergence. Note that denser LVs, when compared to less dense 
ones, can be considered as the result of adding a large-scale varying 
field to the latter. Consequently, we expect that the CW elements 
appear and percolate earlier on within denser LVs than within less 
dense ones. 

These considerations apply to the evolution of the Ilj eigen¬ 
vectors, ei{z), and to their possible dependence on mass. 

Regarding shape evolution, as already emphasised, mass 
anisotropically flows towards new singularities. These anisotropic 
mass arrangements make the eigenvalues evolve. Thus, evolu¬ 
tion becomes gradually extinct as anisotropic flows tend to vanish. 
At small scales, the CW structure is swallowed up and removed by 
contractive deformations, see previous subsection. From a global 
point of view, the CW dynamic evolution somehow stops and the 
structure becomes frozen as ^ 

dominates the expansion at za, see equation]^ Therefore, matter 
flows are expected to become on average less and less relevant af¬ 
ter za, as time elapses. 

In addition, it is expected that locally the first to vanish are the 
flows associated with «(</), the largest eigenvalue of the local de¬ 
formation matrix dij {q} (i.e. the flows towards walls), and the last 
to disappear are those flows related to 7 (g), the smallest deforma¬ 
tion matrix dij{q) eigenvalue (i.e the flows towards nodes). 

Disentangling how these theoretical local predictions affect 
the global shape evolution of LVs demands numerical simulations. 
We will address these issues in the next sections. 


4 EVOLUTION BEYOND THE ZA OR THE AM 

Some concepts, not directly described by the ZA or the AM, need to 
be clarified in order to correctly explain Fig.|^at a qualitative level, 
as well as some results to be discussed in forthcoming sections. 

4.1 Caustic dressing 

The phenomenological Adhesion Model tells nothing about the in¬ 
ternal density or velocity structure of locations where mass gets 
adhered. Just to have a clue from theory, we recall that in his 
derivation of a generalised adhesion-like model, [Domfngue^j2000^ 
found corrections to the momentum equation of the ZA that regu¬ 
larise (i.e., dress) its wall singularities. These then become long- 
lasting structures where more mass gets stuck, but within non-zero 
volumes supported by velocity dispersion coming from the energy 
transfer from ordered to disordered motions, (see also |Gurbatov,| 
jSaichev & Shandarin|1989| for a discussion of these effects in terms 
of the viscosity, phenomenologically introduced in the AM). The 
analyses of A^-body simulations strongly suggest that any kind of 


flow singularity gets dressed (i.e., not only at pancakes, as it has 
been analytically proven by|Dommguez|2000^. 


4.2 Gas in the cosmic web 

When gas is added, the energy transfer from ordered to disordered 
motions around singular structures includes the transformation of 
velocity dispersion into internal gas energy (heating) and pres¬ 
sure. Then, energy is lost through gas cooling, mainly at the dens¬ 
est pieces of the CW, making them even denser. However, as al¬ 
ready said in §3.3| singular (i.e., dense) and regular (i.e., not yet 
involved in singularities, low density) mass elements are mixed at 
any scale. Therefore, low-density gas is heated too, and, in addi¬ 
tion, pressurised. The consequences of these processes cannot be 
deciphered from theory, but previous analyses of cosmological hy- 
drodynamical simulations in terms of the CW (see, for example 
jPomlnguez-Tenreiro et al.|2011^ suggest that dressing acts on any 
kind of flow singularity, i.e., also on filaments and nodes. Moreover, 
these authors conclude that, at (node-like) halo collapse, cooling of 
low-density gas is so slow that most gravitationally heated gas is 
kept hot until 2 ; = 0. In any case, because hot gas is pressurised, 
no anisotropic mass inflows towards singularities can be expected 
within the hot gas component, on the contrary, possible anisotropic, 
pressure-induced hot gas outflows are expected from them. These 
expectations will be explored in the following sections. 

On the other hand, at the densest gas locations, cold gas is 
transformed into stars with an efficiency e when the density is 
higher than a threshold. In this way, the hot gas component and 
the stars, observed in Fig.|^ arise. 


4.3 A visual impression of LV evolution 

Fig. 0 gives us a first visual impression of the evolution of the 
initially spherical LVs. The former considerations above make it 
easier a qualitative interpretation of what these figures show. In¬ 
deed, the gradual emergence of a local skeleton stands out in both 
of them, including web-element mergings and some rotations too. 
Finally, at ziow, we see an elongated structure, either in the DM, 
cold or hot baryonic components, where different spherical con¬ 
figurations appear, with a stellar component at the centre of most 
of thenj^ A high fraction of hot gas component (but not its whole 
mass) is related to these spheres. This complicated structure comes 
from wall and filament formation, according to the AM, and its 
dressing and eventual fragmentation into clumps. Clumps are in 
their turn dressed. Note also that, at each 2 , a fraction of the matter 
is not yet involved into singularities. Therefore, evolution leads to: 
(i) a DM component sharing both a diffuse and a dressed singularity 
configuration, with the particularity that the LV diffuse component 
present at redshift 2 has not yet been involved in any singularity 
at 2 , (ii) a complex cold gas component, sharing also a diffuse as 
well as a dressed singularity configuration, but with a more concen¬ 
trated distribution than that of the DM, because gas can lose energy 
by radiation and (iii) a complex hot gas distribution. As explained 
in § |4.2| diffuse gas is gravitationally heated at collapse events, but, 
as will be shown in § |6.3| it is not involved in important anisotropic 
mass rearrangements. 

To further advance, we need a quantitative analysis of LV evo¬ 
lution. This is the subject of the next sections. 


® We note that there is a component effect, namely different components 
(i.e. DM, cold and hot baryons) evolve dissimilarly. 
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5 ANISOTROPIC EVOLUTION: EIGENVECTORS OE 
THE MASS DISTRIBUTION 

According to the AM, mass elements are anisotropically deformed 
and a fraction of them pass through one or several singularities in 
sticking regions. For each mass element placed at a Lagrangian 
point q, accretion at high jz preferentially occurs along the eigen¬ 
vector corresponding to the largest eigenvalue of the symmetric de¬ 
formation matrix at q, dij (q) = — • 

Taking the LV as a whole, the ir eigenvector 63 °* (z) which 
corresponds to its larger eigenvalue, A 3 ( 2 ) at a given redshift 
z, defines the direction along which the overall LV elongation 
has been maximum until this 2 . Similarly, ei°*( 2 ) corresponds 
to the direction of overall minimum stretching of the LV up to 
a given 2 . It is very interesting to analyse whether or not there 
exists a change in such directions as cosmic evolution proceeds. 
In Fig. 1^ we show the histograms for the quantities Ai{z), the 
angle formed by the eigenvectors and ef^lzio-w), with 

i = 1 , 2 ,3, where ‘tot’ stands for the eigenvectors of the jr ten¬ 
sor corresponding to the total mass of the LV, at redshifts 2 = 
10, 5, 3,1, 0.7,0.5,0.25, 0.25, 0.1. That is, we measure the devi¬ 
ations from the eigendirections at a given 2 with respect to the 
final oneQ We see that on average these directions are frozen at 
Zfvoz ~ 0.5, in such a way that only a few LVs change the eigenvec¬ 
tors of their total mass distribution at 2 ^ Zf^oz, while at 2 ^ 2 froz 
more and more LVs do it. This behaviour is illustrated by Fig.|^ 
where the evolution of the Ai{t) for a typical LV case is plot¬ 
ted. We observe that Ai{z) smoothly and gradually vanish before 
t/t\j = 1, this behaviour being common to all the LVs. 

This is particularly interesting, because as we will see in Figs 
[^and|^the evolution of the jr eigenvalues (or, equivalently, that of 
its principal axes of inertia a, b, c), also declines before t/t\j = 1 . 

It is also important to investigate if there exists a component 
effect in the freezing-out of the eigendirections. With this purpose, 
we have compared the directions of the principal axes of inertia 
that arise from the whole mass distribution with the ones derived 
from every component at different redshifts (see Fig.|^. We have 
found that the latter are mainly parallel to e*°* in the DM and cold 
baryon cases. Concerning hot gas, the distribution of the angles, 
Oi, formed by 6 *°* and the eigenvectors of the hot gaseous 

component, starts nearly uniform and as time elapses a peak around 
0° arises, as we can observe in Fig.|^for the ei case. 

This means that DM dynamical evolution determines the pre¬ 
ferred directions of LV stretching, and cold gas particles closely 
follow them. Hot gas particles (in this case, as explained in Q 
gaseous particles not trapped into singularities and heated by grav¬ 
itational collapse), on the contrary, do not follow DM evolution at 
high redshifts, but they trace at any 2 the locations where mass 
sticking events have taken place. Indeed, as explained in SQ gas 
gravitational heating is due to the transformation of the ordered 
flow energy into internal energy at CW element formation. 


6 ANISOTROPIC EVOLUTION: SHAPES 

Before we focus on the statistical analysis of our results, we present 
the shape evolution of some selected LVs in order to show how 
they acquire their filamentary or wall shape. Then, we analyse the 

^ Note that only two out of the three Ai angles are independent in such a 
way that if for instance Ai = 0 then A 2 = A 3 . 



0 20 40 60 80 20 40 60 80 20 40 60 80 

A, A 2 Ag 


Figure 3. Evolution across redshifts of the Ai distribution, where Ai is the 
angle formed by the eigenvectors el°* ( 2 ) and el°*^( 2 io,„), with i = 1,2,3, 
and where ‘tot’ stands for the eigenvectors of the /E, calculated with all 
the LV components. 


shape evolution of all the objects in our sample, by considering 
component as well as mass effects. To that end, LVs are grouped 
according to their mass, M, into three bins, massive (M ^ 5 x 
IO^^Mq), intermediate mass (5 x 10^^ < M < 5 x 10^^Mq) and 
low-mass LVs (M < 5 x 10^^Mq). 
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z 



Figure 4. An example of the Ai (f) evolution, where Ai is the angle formed 
by the eigenvectors eY'^(z) and with i = 1, 2, 3 and t is given 

in terms of the age of the Universe (fu). 


DM Cold bar. Hot gas 



Figure 5. Distributions of the angles formed, at several redshifts, by the 
direction of the axis of inertia that arise from the overall matter 

distribution with the same axis calculated with the different components. 

6.1 Two particular examples of shape evolution 

In Fig. 1^ we exemplify the evolution of the principal axes of the 
inertia ellipsoid for the LVs of Fig. The upper plot (LV on the 
left-hand side of Fig. illustrates an LV that has two axes that 
expand across time, i.e., it has a flat structure. The lower plot cor¬ 
responds to the LV on the right-hand side of Figure]^ and portrays 



0.0 0.2 0.4 0.6 0.8 1.0 

t/tu 


Figure 6. Evolution of the principal axes of inertia for two LVs. Top, LV on 
the left-hand side of Fig. with a wall-like structure. Bottom, LV on the 
right-hand side of Fig.which acquires a filamentary shape. 



Figure 7. Axis ratio evolution of the Lagrangian volumes of Fig. The 
upper plot shows the evolution towards an oblate shape and the lower plot 
shows an LV that acquires a prolate shape. 

the case in which the major axis grows while the other two axes are 
compressed, giving in consequence a prolate shape. This result can 
also be inferred from Fig. [7] where we can see the evolution of the 
axis ratios h/a and c/a for the same LVs of Fig.In the lower plot 
of Fig.|^ we observe that the two minor axes end up close to each 
other in length, therefore the LV has a filamentary structure. The 
upper plot, in contrast, has the minor axis significantly shorter than 
the other two, hence having an oblate shape. 

A remarkable result is the continuity of the o(f), b{t) and c{t) 
functions for all the LVs, with no mutual exchange of their respec¬ 
tive eigendirections across evolution, i.e., the local skeleton is con¬ 
tinuously built up, in consistency with [Kidding, Shandarin & van| 
|de Weygaert| ( |2014T l. 

6.2 Generic trends of shape evolution 

In this subsection, the generic trends of shape evolution are exam¬ 
ined at a qualitative level. In Fig.[^ where the axis ratios are plot- 
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ted, we can note that the selected LVs are gathered on the nearly 
spherical zone ic/a ^ 0.8) by construction, except the hot gaseous 
component. As time elapses, LVs are deformed, and their evolution 
is shown as they move down inside the triangle described by the 
axes h/a, c/a and T = 1 (orange line). Accordingly, at 2 = 0.05 
they tend to be spread over the triangle. Note that intermediate mass 
and low-mass objects evolve faster than the massive ones. At 2iow, 
DM is preferentially located in the T > 0.3 and c/a < 0.4 region, 
therefore we end up with more prolate systems than oblate objects. 
This assertion is valid for the total, DM and cold baryons axis ratio 
evolution. In contrast, hot gas does not seem to show a remarkable 
evolution effect as it appears populating roughly the same regions 
of the aforementioned triangle at redshifts 10,5 and 3, and later on, 
excluding either the oblate area on the right or the prolate one at 
the left bottom comer of the triangle. 

The shape evolution of the LV mass distribution is also shown 
in Fig.|^ where shape distortions are represented in the prolateness- 
ellipticity plane. In this case, LVs move inside the triangle bound 
by the lines, e = p (prolate spheroids), p = —e (oblate spheroids) 
and p = 3e — 1 (flat objects). We observe the same pattern as in 
Fig. 0 for the total components, DM and cold baryons. In other 
words, initially spherical systems, concentrated on one corner of 
the triangle, evolve across redshifts filling up the triangle, so that, 
at 2 = 0.05, we end up with a high percentage of prolate triaxial 
objects, ~ 83% for the total inertia ellipsoid. We have also found 
that ~ 91% of the selected LVs have extreme total ellipticities (e > 
0.5), while only 8% have moderate ones. A significant percentage 
of the analysed objects are extremely prolate, ~ 31%, that is, they 
have a thin filament-like shape. At 2 = 0.05, we can find systems 
close to the fiat limit, specially in the case of cold baryons. As in 
the previous figure, hot gas does not present a remarkable evolution 
effect after 2 = 1. At higher 2S, however, the hot gas in some LVs 
show needle-like as well as flat shapes (see panels corresponding 
to 2 = 5 and, to a lesser extent, at 2 = 3), but these shapes do not 
appear anymore at lower 2S. 

Figs and nicely show generic trends of shape evolution. 
More elaborated, quantitative analyses of component and mass ef¬ 
fects are given in the next sub-sections. 

6.3 Component effects 

In order to quantitatively determine if there is a component effect 
on the LV shape evolution (i.e., whether DM, hot and cold baryons 
behave dissimilarly), we represent the cumulative distribution func¬ 
tion (CDF) of the e, p and T parameters in Figs.|10|and|l 1[ 

Each row in Fig. [T^ shows the cumulative probability of the 
e parameter calculated for DM, cold baryons, hot gas and the total 
components at a given redshift. The first column depicts the result 
obtained for all the LVs and the other columns display our find¬ 
ings split according to the binning in LV mass. As we can observe, 
the DM and cold baryonic components move from low ellipticities 
or high sphericities at high redshifts towards higher ellipticities at 
2 iow. As a result, these components acquire a filament-like struc¬ 
ture (see Fig.|10^. Note that cold baryons and DM exhibit approx¬ 
imately the same behaviour as time elapses. At 2iow cold baryons 
are slightly more prolate than the DM component, specially in the 
case of low-mass LVs. On the other hand, the hot gaseous com¬ 
ponent does almost not experience an evolution effect, as can be 
noted from the ellipticity CDFs in Fig.[T^ whether or not we group 
the LVs according to their mass. Hot gas has an e ~ 0.57 since 
2 = 2, and does not present any preference for either a spherical or 
a filamentary structure. 


Similar conclusions can be extracted from the DM, cold 
baryon and hot gas prolateness CDFs (see first row of Fig.|l l|l. In 
this case, hot gas has anpranging from 0.25 —0.34 since 2 = 2. An 
important difference with respect to the ellipticity CDFs is that at 
2 iow, hot gas cumulative probabilities show a small deviation from 
cold baryons CDFs which is bigger in the low-mass bin, while in 
the e case these components exhibit a large deviation from each 
other. 

Triaxiality CDFs show a tendency of cold baryons to have a 
prolate shape independently of the mass binning at 2 = 3. We ob¬ 
serve the same displacement of DM and cold baryon CDFs across 
redshifts, previously noted from ellipticity and prolateness cumu¬ 
lative probabilities. Concerning hot gas, it has an T in the range 
0.69 — 0.76 since 2 = 2, showing almost no changes thereafter. 
This displacement causes that the difference between DM, cold and 
hot baryons CDFs appears greatly diminished at 2 = 1. This fact 
can also be noticed from ellipticity and prolateness CDFs. It is note¬ 
worthy that the cold baryon triaxiality cumulative probability of the 
massive LV bin is delayed with respect to the DM CDF at 2 = 1. 
This difference is kept at 2 = 0.05 (see lower panels of Fig. \n}^ 
this is also true for the prolateness case. On the contrary, at 2iow 
DM CDF appears delayed with respect to cold baryons for the low- 
mass bin. 

6.4 Mass effects 

To study the impact of the LV mass on its shape deformation, we 
plot in Figs[T^and [T^ the CDF split by the component considered 
in the reduced inertia tensor calculation. From left to right we show 
in each column results obtained with all the particles, taking into 
account only DM particles, then cold baryons results and finally 
hot gas. Rows in Fig.[T^show cumulative probabilities at different 
redshifts. Each panel present the CDF calculated according to the 
binning in LV mass, massive object CDF are shown in magenta, 
intermediate mass results in cyan and low-mass CDF in blue. 

In the first place, we discuss the ellipticity CDFs in Fig.[T^ As 
we can observe, the mass effects are not very relevant and moreover 
they almost do not evolve. The most important mass effects appear 
in cold baryons at any 2 . Indeed, the massive and low-massive LV 
samples at 2 = 3 and 1 have been determined to be drawn from dif¬ 
ferent populations with the two-sample Kolmogorov-Smirnov test 
with 90% Cl; while the massive and intermediate mass LV samples 
with 95% Cl at 2 = 3,1 and 0.05. In general, massive LVs tend 
to be more spherical across redshifts, and they have a narrower e 
distribution than less massive ones. 

In the prolateness case, the mass effects grow with time, ex¬ 
cept in the hot gaseous component. Hot gas independently on the 
mass binning is less spherical than the other components at 2 = 3. 
At 2 = 1, massive LVs are more spherical than the less massive 
ones for both DM and cold baryons. The mass effect is less pro¬ 
nounced in the case of hot gas. At 2iow the tendency described 
above is kept (see upper panels in Fig. |13|l. The p distribution in 
massive LVs is narrower than those in the other mass bins and it 
becomes wider faster in the low-mass bin. 

Regarding triaxiality CDFs, again mass effects grow with evo¬ 
lution, mainly in the DM component (see lower panels in Fig.|13|l. 
We can also note that in both, the total and the DM case, there 
are almost no systems with T < 0.6, specifically, there is a lack 
of oblate massive objects relative to the other mass groups. We 
have tested the difference between the massive and the low-massive 
bins with the two-sample Kolmogorov-Smirnov test at a 90% CL 
This mass effect is less significant in the baryon case. Indeed, cold 
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b/a b/a b/a b/a 

Figure 8. Axis ratio evolution of all the selected LVs, where coloured circles indicate different mass range. Massive LVs with M ^ 5 x IO^^Mq are 
represented in red, LVs with intermediate mass, 5 x 10^^ ^ M < 5 x in cyan and low-mass LVs, M < 5 x in blue. The orange line 

con'espond to T = 1, i.e., to a prolate spheroidal shape. Objects with c/a < 0.9 and T > 0.7 (magenta line) have a prolate triaxial shape and LVs with 
c/a < 0.9 and T < 0.3 (green line) are prolate triaxial ellipsoids. We show the axis ratios obtained with the total number of particles, the axis ratios of DM 
particles, and the axis ratios found for cold and hot baryons. 
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Figure 9. Prolateness-ellipticity plane for the reduced inertia tensor of the selected LVs for redshifts 10,5,3,1 and 0.05. Massive LVs with M ^ 5x10^27^^^ 
are represented in red, LVs with intermediate mass, 5 x 10^^ ^ M < 5 x IO^^Mq, in cyan and low-mass LVs, M < 5 x IO^^Mq, in blue. The orange 
lines correspond to ultimate shapes, e = p (prolate spheroids), p = —e (oblate spheroids) and p = 3e — 1 (flat objects). 
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Figure 10. Cumulative distribution function of the ellipticity parameter, e, portraying component effects and their evolution in different mass bins. Each 
column shows the distribution binned according to the LV mass. Plots in the first column are calculated for the total number of LVs. Rows represent different 
redshifts. The code colour used in each plot is as follows, results obtained with the total reduced inertia tensor are presented in blue, DM results in magenta, 
cold baryons in green and hot gas in red. 




Figure 11. Upper panels, CDF of the prolateness parameter, p, at ziow Lower panels, CDF of the triaxiality parameter, T, parameter at 2 : 10 ^. Each column 
shows the distribution binned according to the LV mass. Plots in the first column are calculated for the total number of LVs. The code colour is as in Pig.|lQ| 


baryons do not present a significant mass effect, only less massive 
LVs tend to be more oblate than the more massive bins at 2 = 3 
and 1. 

Summing up, except for the hot gas component, more massive 
LVs tend to evolve slightly more slowly from their initial spherical 
shape than less massive ones. This can be interpreted in terms of 
the CW dynamics as follows: more massive objects would appear 
more frequently in nodes of the CW, versus less massive objects 


being present in filaments and walls. Therefore, the relative im¬ 
portance of anisotropic mass rearrangements versus radial ones is 
lower in massive than in less massive LVs. Concerning the hot gas 
component, no relevant evolution has been detected, particularly 
after 2 ~ 3, indicating that neither the possible anisotropic flows 
towards singularities, nor the possible pressure-induced anisotropic 
outflows, have caused measurable LV mass rearrangements in the 
LV sample thereafter. 



















Lagrangian Volume Deformations 


15 



e e e e 


Figure 12. Cumulative distribution function of the ellipticity parameter, e, illustrating mass effects and their evolution according to the LV components. Each 
column displays the distribution binned according to the components taken into account to calculate the reduced inertia tensor, namely, the total number of 
particles, DM, cold baryons and hot gas. Rows represent different redshifts. In each plot, massive LVs (M 5 X IO^^Mq ) are shown in magenta, LVs with 
an intermediate mass (5 X 10^^ Si M < 5 X 10^^Mq) in cyan and low-mass LVs (M < 5 X in blue. 




T T T T 


Figure 13. Upper panels, CDF of the prolateness parameter, p. Lower panels, CDF of the triaxiality parameter, T. From left to right the columns show the 
distribution binned according to the components taken into account to calculate the reduced inertia tensor, i.e., the total number of particles, DM, cold baryons 
and hot gas. The code colour in each plot is as in Fig.|12| 


7 FREEZING-OUT OE EIGENDIRECTIONS AND 
SHAPES 

7.1 Ereezing-out times 

In the previous sections we have become aware that the Ai{z),i = 
1, 2 and 3 angles evolve with time and —> 0° before ^low We re¬ 
mind that Ai{z) is the angle formed by the eigenvectors el°*{z) 


and ei°*(2iow), with i = 1,2, 3, where ‘tot’ stands for the eigen¬ 
vectors of the lij tensor corresponding to the total mass of the LV. 
Also, the evolution of the LV inertia ellipsoid declines in the same 
limit, see Figs |^and In this section, we use the times when 
these eigendirections and inertia axes become frozen. We have cal¬ 
culated these freezing times to study and compare both processes 
and to look for possible mass effects. The subject is interesting to 
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Figure 14. Histograms for and t™’" defined with 

cos(SAi) = 0.9 and / = 0.1. 



Figure IS. CDFs for the same quantities in the previous figure, showing 
possible mass effects. 


elucidate how and when the local CW around galaxies-to-be be¬ 
comes frozen at the scales analysed in this paper, while it still feeds 
the protogalaxies at smaller scales. 

Having the Ai(z) angles ~ 0° during a « range 2 ^ ziow 
means that the LV deformations become fixed in their eigendi- 
rections before 2iow, or, in other words, mass rearrangements are 
thereafter organised in terms of frozen symmetry axes making the 
inertia tensor diagonal, i.e., in terms of a skeleton-like structure. 
This motivates the search for the moment when a given LV gets its 
structure frozen. This is not a straightforward issue, however, be¬ 
cause this situation is gradually reached: all we can do is to resort 
to thresholds. 

In the following, we use time instead of 2 in order to make 
our results clearer. Given a threshold angle SAi, we define tsAi 
as the time (Universe age at the event in units of the current Uni¬ 
verse age tu) when Ai(t) ^ SAi if t ^ tsA,, (i.e., the Universe 
age when the ith eigendirection of the inertia tensor becomes fixed 
within an angle SAi). Then, we define as the maxi¬ 

mum and minimum values of tsAi ,i = 1,2,3, for each LV. That is, 
t^A^ for a given LV is the fractional time when the directions of its 
three eigen vectors become frozen, or, symbolically, Ai(t) ^ SAi 
if t ^ t^A^ for any directioij^ The minimum satisfies the 
same condition for just one direction. Fig. |14| (upper plots) shows 
the distribution of and for our sample of 206 LVs with 
SAi such that cos((5Ai) = 0.9. 

A very interesting point is to explore LV shape transforma¬ 
tions relative to the freeze-out times for inertia eigendirections. An 
illustration can be found in Figs and Comparing both figures, 
we see that the principal axes change slightly after skeleton emer¬ 
gence for the particular LVs considered in this figure by using a 
10% threshold (see below). The differences are larger for other 
LVs, and, indeed it is worth analysing this issue in more detail. 


® Note that the second and the third eigendirections become frozen at the 
same time. 


Therefore, to be more quantitative, we define as the fractional 
time when the inertia axis a becomes frozen within a threshold 
fa, which is a fixed fraction of the a{t) value, i.e. Aa{t) fa 
if t ^ tf^a, where Aa{t) = ■ Similarly, we define 

tf.f, and tf^c, and then and The former is the time when 
the three inertia axes become frozen, while the latter is the time 
when just one axis gets frozeij^ To have an insight of the statistical 
behaviour of these times, in Fig. (lower plots) the histograms 
for and are represented for / = 0.1. In this figure, right- 
(left-)panels correspond to the times when one (three) out of the 
eigenvectors or the principal inertia axes become fixed within a 
10 % of their final values. 

An interesting result is that the time range for is narrow 
and late. The range of t^X^ ** much wider, which means that a 
high fraction of LVs get at high 2 their three eigendirections fixed 
before the evolution of their inertia axes ends up. During this early 
time interval, LVs change their shape with frozen symmetry axes, 
i.e., anisotropic matter inflows onto CW elements. Another result 
is the accumulation at the first bin of the evolution time: these 
are the systems having a principal axis of inertia that keeps within 
a 10% of its initial value along the evolution. They are less prolate 
than other systems. An even higher fraction of LVs have one of 
their eigendirections fixed in the first 5% of the evolution time (see 

Fig.frgb). 

A high fraction of systems also got one frozen eigendirection, 
while none of their principal inertia axes is fixed yet. However, at 
the end of the evolution this effect vanishes (compare Figs |14| b 
and|14[d). Finally, let us mention that LVs also spend an important 
fraction of their lives with one but not three fixed eigendirections 
(within the thresholds used to draw these figures, compare Figs|14|a 
and|14|b), or one but not three frozen inertia axes (compare Figs 
[T4 ]c and[T4]d). 

® Again, once the value of one principal axis becomes fixed, the freezing 
times for the other two axes are the same. 
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7.2 Mass effects 


Next, we look for mass effects in the distributions of and 
as well as in those of and f7““. This is more clearly visualised 
in terms of cumulative histograms. In Fig.[^ we plot the CDF for 
and (i.e., LV eigen directions relative to their final values, 
first row) and and (principal inertia axes, second row), 
respectively, where no binning has been used. To analyse possible 
mass effects, results for the three mass groups are shown in each 
panel. The cumulative histograms in the four panels of this figure 
are in one-to-one correspondence with the histograms in Fig.|14[ 
The first outstanding result is that the time range for 
is roughly the same (narrow and late), irrespectively of the mass 
range used (Fig j^l^|c). This behaviour can be understood as the 
consequence of ——>■ 0 at late times, a global effect causing 
anisotropic flows to vanish, see §3.4| for more details. Nevertheless, 
there exists a mass effect in (Fig.[^a), with the least-massive 
LVs showing a delay in the spine emergence or in getting their three 
eigendirections frozen with respect to more massive ones, the dif¬ 
ferences being more marked at early times. This is somewhat ex¬ 
pected from the previous discussion on the effects of the eigenvalue 
landscape heights on the timing of spine emergence, in p.4| 

Fig. [T^b exhibits strong mass effects too. Indeed, at early 
times the most massive systems get one out of their three eigendi¬ 
rections frozen sooner than less massive ones. In fact, ~ 95% of 
the massive LV subsample has one of their eigendirections fixed at 
t/t\j ~ 0.1. This mass segregation can be understood in the light 
of the considerations made in §3.4[ where we concluded that the 
first CW elements tend to appear and percolate earlier on within 
massive LVs than within less massive ones. 

On the other hand, the freezing-out times for the principal axis 
of inertia (panel [T^d) display a remarkable mass effect, although 
just at early times. Later on, irrespective of their mass, no LV gets 
its first principal axis of inertia fixed later than t/tu ~ 0.55. 
This upper bound on might be a consequence of both, the 
0 after the A term dominates the Universe expansion, 
and the fact that flows towards walls are the first to vanish at a lo¬ 
cal level. The mass effect lies in massive systems having their tf™ 
delayed at early times in relation to less massive ones (consistently 
with what was found in § |6.2^ , the difference vanishing at 2 ~ 1. 

Finally, to look for correlations, the and for our 
sample of LVs are plotted versus their respective 


in Fig. 16 for / = 0.1 and cos(5Ai) = 0.9. No outstanding corre¬ 


lation exists in any case, but we see that indeed, most systems have 
their eigendirections fixed before their principal axes got frozen. 

Summing up, we observe that on average eigendirections (ei¬ 
ther one or the three) for massive LVs become fixed at earlier stages 
than that of less massive LVs. Nevertheless, no relevant mass ef¬ 
fects are found for principal inertia axis freezing times. In addition, 
eigendirections become in general fixed before mass flows onto the 
corresponding CW elements stops, the time delay being particu¬ 
larly long for the first eigendirection relative to the first principal 
axis in massive systems. Thus, the first eigendirection in massive 
systems gets fixed quite a while before the accretion onto it stops. 



Figure 16. Scatter plots of versus (Lft) ^nd t'?"' versus t™’ 
(right). 


of LVs in the sample, and a high K, ensuring LVs with high enough 
number of particles so that we obtain meaningful LVs. However, 
K — 10 is by no means the unique value that satisfies these con¬ 
straints. Therefore, it is important to test out the possible effects of 
changing this value under the same constraints. 

To this aim, we have repeated all the calculation using K = 
7.5 and 15. The LV building up (see section 2.2) has been repeated 
with the same SKID identified haloes at ziow as first step. Nonethe¬ 
less, when iV = 15 is used, some of the LVs do not satisfy anymore 
the condition of having all their particles inside the hydrodynamic 
zoomed volume. These particular LVs have been removed from the 
initial sample of 206 LVs, in such a way that we are finally left 
with 159 LVs for K = 15. This problem does not exist when us¬ 
ing K = 7.5; however, to probe the scale effects, we need samples 
that contain the same ziow SKID-identified haloes as starting point 
in the three scales. Therefore, only these 159 well-behaved LVs (a 
subset of the initial K = 10 sample) have been used to analyse the 
scale effects. 

The first relevant outcome is that there is no substantial dif¬ 
ference when results obtained with the subsample of 159 LVs and 
with the sample used along this paper (206 LVs) for K — 10 are 
compared. 

In the following subsections, we will compare the results ob¬ 
tained with each of the three samples of 159 LVs, dubbed according 
to its K value, ify.s, Kio and TTis. 


8.1 Effects on eigenvector orientation evolution 

Concerning the evolution across redshifts of the jr eigendirections 
relative to their final values at z\o-w (Fig.[^, no relevant differences 
have been found between the histograms obtained with the TTis and 
K\o samples at the same redshifts. Fig.|17|illustrates this behaviour, 
showing that the Ai angle distributions for Tfis are similar to those 
found with Kiq at different 2 pairs, see §8.3| for more details. 

In addition, no scale effects appear in the angles formed by the 
eigenvectors, i = 1,2,3, arising from the overall matter 

distribution with the same eigenvectors calculated with the different 
components (i.e., those angles whose distribution for the sample of 
206 objects is given in Fig.|^) 


8 DISCUSSION: POSSIBLE SCALE EFEECTS 

In Section [2^ when describing how to build up the LV sample, a 
value of iihigh = K X Tvir.iow with K = 10 has been chosen to de¬ 
fine the LV at 2high. As explained there, this choice was motivated 
as a compromise between low K values, ensuring a higher number 


8.2 Effects on shape evolution 

To gain further insight, the 159 LV subsample has been split ac¬ 
cording to the LV masses. In order to assure that we are comparing 
the same mass bins for the three scales, we have mapped the LVs 
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Figure 17. Histograms of the Ai distribution at different redshifts for the 
Kio and K 15 samples (left- and right-hand columns, respectively). 


belonging to the three mass ranges defined for the Kio sample to 
the LVs of the Ki^ and Jfr.s scales. 

Important results concerning shape evolution are as follows. 

(i) No relevant differences in the evolution patterns have been 

found between the least massive LV group (M < 5 x in 

the Kio sample) when followed in the K 15 , Kio and ifr.s samples 
(see Fig. [T^ blue lines). That is, these LVs are hardly sensitive 
to the K scale in their evolution. The scale effects are only slight 
between the K 15 and Kiq samples when no mass splitting in the 
LV sample is performed (see Fig.[T^ black lines). 

(ii) LVs in the massive group are sensitive to the K scale, with 
the iTr s samples showing particular differences. Fig.[T^is an ex¬ 
ample of such a behaviour, likely due to the wall effect, whose for¬ 
mation is better sampled with Jfis. Also, walls are more frequent 
in massive LVs. See j ]8.3| for more details. 

(hi) In any case, the qualitative results reached in §|^about com¬ 
ponent effects in shape deformations are stable when comparing 
fTis and Kio samples. 
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Figure 18. CDFs of the ellipticity at portraying mass effects obtained 
with the three different scales, Ar.s Kio and K 15 . 



t“'”(tu) tp'”(tu) 

Figure 19. Histograms for and t™*" defined with 

cos(SA-i) = 0.9 and / = 0.1. Columns show the results obtained for 
the three samples, Ky.s Kio and Kis- 


8.3 Effects on freezing-out times 
Fig, [r?! shows the histograms for the 

times for samples using different K scales. It is clear from this fig¬ 
ure that while the results for the Kiz and Kio samples are roughly 
consistent with each other, those for the Kj,^ sample differ. The 
only exception is the time distribution (second row), whose 
pattern is the same at any scale, namely rather late and peaked. Re¬ 
call that is the time when the three inertia axes are fixed to 
within 10% of their final values, i.e., the time when all anisotropic 


fluxes stop. This behaviour can be understood as the consequence 
of > 0 at late times, that is a global effect. 

A key point to understand some aspects of Fig.|19|behaviour, 
is the fact that the iTy s scale is too short to suitably sample the 
whole process of wall formation within some LVs. As a conse- 
quence, since the first flows to vanish are those towards walls (see i 


3.4 1 , the f™'" time (when the first inertia axis is fixed to within 10% 


of its final value) will be delayed at high 2 in the Kr .5 sample, as 
observed in Fig. fourth row. A remarkable results is that, irre¬ 
spective of the K scale, no LV has its first inertia axis frozen later 
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Figure 20. CDFs of tff, and at different K scales. 


than t/t\j ~ 0.55. This result reinforces our interpretation given 
in § |7.2|that this effect is, at least partially, a consequence of the 


0 tendency at latter times. 

The process of wall formation could be also the reason of the 
similarities and differences found in the distributions of the 
times (when the first eigenvector direction is fixed to within a 10 %). 
The panels of the third row of Fig.[T^show that their distributions 
are always peaked towards very early times, meaning that the 63 
eigenvector of the ir for some LVs freezes its direction very early, 
following wall formation. In addition, we see that as we move from 
Tfis to Kit) to Jfr.s, a delay appears, not so relevant between the 
fTis and Kiq samples. Again, this can be interpreted in terms of the 
inadequacy of the shorter scale to properly catch the characteristics 
of wall formation in some LVs. 

Finally, we address the scale effects on (first row of 
Fig. [T^ . These are the times when the LV orientations become 
frozen to within 10% of their final values, i.e., the times marking 
the skeleton emergence locally within each LV. While its distribu¬ 
tion is rather peaked at early times for both, the K 15 and the Kio 
samples, it flattens as we go to Kta- Once again, the poor wall for¬ 
mation sampling in most K 7.5 LVs is likely to be the cause of this 
difference. 

It is worth noting that the qualitative features found in §|^are 
stable under the change in K. For instance, mass effects can be 


analysed from Fig. 20 where we show the and 


f™'" mass-binned CDFs (first, second, third and fourth rows, re¬ 
spectively), at different scales (columns). Then, we can note that, 
regardless of the K value, the t^A^ and t^A^ distributions show 
qualitatively similar mass effects, with the most massive LV group 
fixing either one or their three eigenvalues earlier on than LVs in 
the intermediate or less massive group (as expected). Moreover, 
the distribution does not show relevant mass effects whatever 
the considered scale. Finally, irrespective of the scale, the dis¬ 
tributions do not show relevant mass effects after t/tu — 0.4, as 
expected from the previous analyses. At low z, some mass segrega¬ 
tion is found, and furthermore, it qualitatively depends on the scale. 
This is the only one exception to the stability under the change in 
K. These results could reflect the difficulty of catching the end of 
the mass flows in only one direction when the contribution of wall 
formation is combined with mass effects. 

Summing up, the differences in the freezing-out times are not 
very relevant when using the K 15 or the Kio samples. Their distri¬ 
butions show similar patterns, in particular when mass effects are 
considered. 


9 SUMMARY AND CONCLUSIONS 


In this paper, we present a detailed analysis of the local evolution 
of 206 Lagrangian Volumes (LVs) selected at high redshift around 
proto-galaxies. These galaxies have been identified at z\avr = 0.05 
in a large-volume hydrodynamical simulation run in a ACDM cos¬ 
mological context and they have a mass range 1 —1500 x 
We follow the dynamical evolution of the density field inside these 
initially spherical LVs from ^high = 10 up to = 0.05, wit¬ 

nessing mass rearrangements within them, leading to the emer¬ 
gence of a highly anisotropic, complex, hierarchical organisation, 
i.e., the local cosmic web (CW). Indeed, at ziovj LVs acquire over¬ 
all anisotropic shapes as a consequence of mass inflows onto sin¬ 
gularities along cosmic evolution, in such a way that some relevant 
aspects of these mass arrangements can be described in terms of 
the reduced inertia tensor Ilj evolution, as given by its principal 
directions and inertia axes, a ^ 6 ^ c. 

Our analysis focuses on the evolution of the principal axes of 
inertia and their corresponding eigendirections, paying particular 
attention to the times when the evolution of these two structural 
elements declines. In addition, mass and component effects (either 
DM, cold or hot baryons) along this process have also been inves¬ 
tigated. 

In broad terms, we have found that local LV evolution follows 
the predictions of the Zeldovich Approximation (ZA, ZeTdovicTil 
1970[ l and the Adhesion Model (AM, [Gurbatov & Saichev|1984t 


Gurbatov, Saichev & Shandarin|[T989| [Shandarin & Zeldovich] 


1989[ Gurbatov, Malakhov & Saich ev|1992[| Vergassola et al. 11994] ) 

when both caustic dressing ( |Dommguez| 2000 1 and mutual gas ver¬ 

sus CW effects (see Section 3 and Domlng uez-Tenreiro et al.|201 
[Metuki et al.|2015) are taken into account. Evolution also entails 
baryon transformation into stars inside the densest regions of the 
web and gravitational gas heating following the collapse. More 
specifically, these are our main results. 

Dark matter dominates dynamically the LV shape deforma¬ 
tions over the baryonic component, as expected from hierarchi¬ 
cal structure formation. Deformations transform most of the ini¬ 
tially spherical LVs into prolate shapes, i.e. filamentary structures, 
in good agreement with previous findings i Aragon-Calvo, van de| 


|Weygaert & Jones|2010[|Cautun et al.|2014^ . Cold baryons follow 

DM behaviour in general, but with some departures from it, de- 
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partures that rise as evolution proceeds. Accordingly, the number 
of LVs having their cold haryonic principal axes in directions that 
differ from the ones calculated with their DM content is negligible 
at Zhigh, and it keeps low along the evolution, but increases with 
time (~ 25% at ziow). On the contrary, the hot gas eigendirections 
have a flatter distribution at Zhigh and then they tend to converge 
to those calculated with DM. However, only ~ half of them reach 
such convergence at aiow This tendency towards convergence is 
due to the fact that the hot gaseous component traces the locations 
where sticking events, in particular filament and node formation, 
have taken place. The mass fraction involved in these processes in¬ 
creases with evolution, and consequently we expect a tendency of 
the hot gas to be aligned with the total eigendirections. 

In terms of shape evolution, a clear component effect has 
been found regarding the way how the evolution occurs. In fact, 
hot gas shapes do not exhibit important evolution because, as said 
above, gravitationally heated gas marks out the places where stick¬ 
ing events have taken place, and because, in addition, no evidence 
for important anisotropic mass rearrangements in this component 
have been found in this paper. The only remarkable effect is that 
the needle-like or flat shapes shown by hot gas in some LVs around 
z — 5, are transformed at lower zs. As mentioned before, DM 
and cold baryons shapes do evolve, with cold baryons achieving 
an even more pronounced filamentary shucture than DM ones as 
a consequence of dissipation. Additionally, some mass effects have 
also been found in the generic evolution of shapes, with lower mass 
LVs evolving towards more pronounced filamentary structures on 
average and earlier on than the more massive ones. 

A remarkable result of our analyses is that the evolution of 
LV deformations declines. This means that both the LV eigendirec¬ 
tions, as well as their principal axes of inertia (a, b and c) values 
become roughly constant before Ziow This is a smooth effect that 
can be only defined in terms of thresholds. Taking a 10% of the final 
values, shape (i.e., a, b and c values) freezing-out time distribution 
has a narrow peak (~ 0.2 at each side) around t/t\j = 0.8. This 
happens later than the freezing-out times for the three LV eigendi¬ 
rections, whose distribution peaks around t/tu = 0.1 and then it is 
flat until t/t\j ~ 0.8 when it decays. 

By plotting individual freezing times for shapes and eigendi¬ 
rections, respectively (see Fig.[T^a), we note that first, most of the 
LVs fix their three axes of symmetry (like a skeleton), and later 
on their shapes are fixed. This result is in good agreement with 


|van Haarlem & van de Weygaert|(|1993|l;|van de Weygaert & Bond 
(|2008t;|Cautm et al.^Mry and |Hidding, Shandarin & van de Wey 


Igaert ( |^14^ findings. Moreover, the ZA and the AM predict that 
walls, filaments and nodes undergo mass flows from underdense 
regions to denser environments, that continue after skeleton emer¬ 
gence. 

As a general consideration, it has been found that mass re¬ 
arrangements at the scales taken into account have always been 
highly anisotropic. Therefore, the mass streaming towards walls 
and filaments has been extremely anisotropic, and, to a lesser ex¬ 
tent, towards nodes as well. In particular, galaxy systems form in 
environments that have a rigid spine at scales of a few Mpc, from 
whose skeleton a high fraction of mass elements that feed proto¬ 
galaxies are collected. 

Due to anisotropic mass accretion, it turns out that in general 
the direction of just one of the LV eigen vectors or the value of one 
of their axes get frozen while the other two still continue changing. 
Again, for each LV there is a time delay between the moment when 
the first of its eigendirections get fixed (happening within the first 
20% of the Universe age) and the moment when the value of one of 


its principal axes becomes constant (peaking around tjtij = 0.35). 
Therefore, we again find a situation where first the flow direction 
is fixed (as a first piece in the skeleton emergence) while the mass 
flows persist. 

Even more interesting because of its possible astrophysical 
implications (see discussion below) is our finding that more mas¬ 
sive LVs fix their skeleton earlier on than less massive ones, either 
considering just one or the three eigendirections. These results are 
not surprising since the dynamical processes involved in the spine 
emergence are faster around massive potential wells. 

Concerning shape transformation decline, there are no rele¬ 
vant mass effects as far as the complete shape freezing-out is con¬ 
sidered. When just one axis value is taken into account, however, 
an early delay of more massive LVs compared to less massive ones 
clearly stands out, delay that vanishes at half of the Universe age. 

When building up the LV sample at Zhigh a value of i?high = 
K X Tvir.iow with K = 10 has been used to define the LV at this 
redshift. This choice was motivated as a compromise between low 
K values, ensuring a higher number of LVs in the sample, and a 
high K, ensuring that LVs are large enough to meaningfully sam¬ 
ple the CW emergence around forming galaxies. As this K = 10 
value is not the unique value satisfying these constraints, the com¬ 
plete analysis has been repeated using K = 7.5 and 15 instead. We 
have found that when using the K = 15 or the K — 10 samples, no 
relevant differences in the LV eigenvector orientations, shape defor¬ 
mations and freezing-out times appear. Therefore, using K = 10 
is in a sense the best choice. 

It is important to remark that no explicit feedback has been 
implemented in the simulations analysed here, but SF regulation 
through the values of the SFR parameters. We remark that the is¬ 
sues discussed in this paper entail considerably larger character¬ 
istic scales than the ones related to stellar feedback. Hence, it is 
unlike that the details of the star formation rate, and those of stellar 
feedback in particular, could substantially alter the conclusions of 
this paper, at least at a qualitative level. Concerning the inner halo 
scale, we recall that to properly explore the impact of SNe feedback 
into filamentary patterns, high enough resolution in order to resolve 
SNe remnants into the Taylor-Sedov phase are needed. Such sim¬ 
ulations are available (the NUT simulations, at sub-parsec scale), 
but only up to 2 = 9 jPowell, Slyz & Devriendt|201 1| |. Therefore, 
we still have to wait to properly understand how SNE feedback 
can possibly affect the CW emergence and dynamics. However, the 
findings so far, at high z, suggest that the filamentary patterns are 
essentially untouched by SNe feedback ( [Powell et al.|2013) . 

9.1 Astrophysical Implications 

The results summarised so far could have important implications in 
our understanding of galaxy mass assembly, raising different inter¬ 
esting issues. 

According to our results, it takes longer for less massive sys¬ 
tems to fix their spine, possibly making it easier for these systems 
to acquire angular momentum through filament transverse motions 
relative to the galaxy haloes. In fact, recent studies on galaxy for¬ 
mation [Kimm et al.|201 rjjPichon et al.|20TT]|Tillson et al.|2015[ 
[Dubois et al.|2014] ( in the CW context, underline the role that fil¬ 
ament motions in the protogalaxy environment could have had in 
endowing filaments, and eventually the adult galaxy, with angu¬ 
lar momentum. If real, this effect could contribute to the mass- 
morphology correlation (see for instance [Kauffmann et al.[2003[ (. 

Our results also point towards (major) mergers events having a 
high probability to occur within filaments. This is an important is- 

























Lagrangian Volume Deformations 21 


sue, though heyond the scope of this paper. In fact, if confirmed, 
this could decrease the allowed merger orbital parameter values 
(see for example, |Lotz et ^|2010[ |Barnes|[201 1^ , as most merg¬ 
ers would have these parameters constrained within the filament. 

Another issue concerns the use of close pairs in merger rate 
calculations from ohservational data, under the hypothesis that 


these systems are bound and about to merge (see, for instance Pat- 


ton et al.|2000 [Bell et al.|2006[ Kartaltepe et al.|2007[ Patton & At- 


field 

2008 

[Robaina et al.|2010 Tasca et al.|2014| Lopez-Sanjuan 

etal. 

2015 

1 . In this respect, some interesting efforts have been made 


to correct the statistics of pairs that are close in angular distance 
from chance superposition effects on the line of sight, (see e.g., 
[Kitzbichler & White|2008| [Patton & Atfield|2008t , whose results 
are used by other authors in this field. Our results reinforce the 
need for these analyses, in the sense that a detailed determination of 
these corrections, including their dependence on the galaxy prop¬ 
erties, merger parameters and environment, could be crucial for a 
more elaborated understanding of the relationship among close pair 
statistics and merger rates. 

Finally we very briefly address the question of the warm-hot 
gas distribution at intermediate scales. Our results point to the web 
structure being marked out by hot gas from high redshifts. Indeed, 
at scales of 4 — 8 Mpc and at Ziow, hot gas traces the CW ele¬ 
ments. Note that there is observational evidence of warm-hot gas 
at large scales in a filament joining Abell clusters A222 and A223 
( [Werner e t al.|20 0^, w here the DM component has also been de- 
tected joietrich et al.[[2012) , and more recently preliminary evi¬ 
dence of hot gas in cluster pairs has been found from the redMaP- 
Per catalogue jRykoff et al.|2014l ( along the sightline of a QSO by 
[Tejos[ ( [2014[ (, (see also his presentation in The Zeldovich Universe, 
Genesis and Growth of the Cosmic Web, 2014, lAU Symposium). 
Our results concern smaller scale structures, and they indicate that 
hot gas traces the CW since the moment when gas is heated at high 
redshift. Indeed, hot gas maps out the sites where the most violent 
dynamical events have occurred, such as filament, and, more par¬ 
ticularly, node formation. Confirming warm-hot gas in filaments at 
different scales is a major challenge for the advance of our under¬ 
standing of galaxy formation (see for example [Kaastra et al.|2013| 
for details). 


ACKNOWLEDGEMENTS 

We thank Arturo Serna for allowing us to use results of simula¬ 
tions. We thankfully acknowledge to D. Vicente and J. Naranjo for 
the assistance and technical expertise provided at the Barcelona Su¬ 
percomputing Centre, as well as the computer resources provided 
by BSC/RES (Spain). We thank DEISA Extreme Computing Initia¬ 
tive (DECI) for the CPU time allowed to GALFOBS project. The 
Centro de Computacion Cientiffca (UAM, Spain) has also provided 
computing facilities. This investigation was partially supported by 
the MICINN and MINECO (Spain) through the grants AYA2009- 
12792-C03-02 and AYA2012-31101 from the PNAyA, as well as by 
the regional Madrid V PRICIT programme through the ASTRO- 
MADRID network (CAM S2009/ESP-1496) and the ‘Supercom- 
putacion y e-Ciencia’ Consolider-Ingenio CSD2007-0050 project. 
SR thanks the MICINN and MINECO (Spain) for financial support 
through an FPU fellowship. 


REFERENCES 

Altay G., Colberg J. M., Croft R. A. C., 2006, MNRAS, 370, 1422 
Aragon-Calvo M. A., Jones B. J. T, van de Weygaert R., van der 
Hulst J. M., 2007a, A&A, 474, 315 
Aragon-Calvo M. A., van de Weygaert R., Araya-Melo P. A., 
Platen E., Szalay A. S., 2010, MNRAS, 404, L89 
Aragon-Calvo M. A., van de Weygaert R., Jones B. J. T., 2010, 
MNRAS, 408, 2163 

Aragon-Calvo M. A., van de Weygaert R., Jones B. J. T., van der 
Hulst J. M., 2007b, ApJ, 655, L5 
Aragon-Calvo M. A., Yang L. E, 2014, MNRAS, 440, L46 
Arnold V. I., 1983, Uspekhi Mat. Nauk, 38, 77 
Bailin J., Steinmetz M., 2005, ApJ, 627, 647 
Bardeen J. M., Bond J. R., Kaiser N., Szalay A. S., 1986, ApJ, 
304, 15 

Barnes J. E., 2011, MNRAS, 413, 2860 

Bell E. R, Phleps S., Somerville R. S., Wolf C., Borch A., Meisen- 
heimer K., 2006, ApJ, 652, 270 
Bond J. R., Kofman L., Pogosyan D., 1996, Nature, 380, 603 
Bond N. A., Strauss M. A., Cen R., 2010a, MNRAS, 406, 1609 
Bond N. A., Strauss M. A., Cen R., 2010b, MNRAS, 409, 156 
Bryan G. L., Norman M. L., 1998, ApJ, 495, 80 
Buchert T., 1989, A&A, 223, 9 
Buchert T., 1992, MNRAS, 254, 729 
Buchert T., Dommguez A., 1998, A&A, 335, 395 
Buchert T., Dommguez A., 2005, A&A, 438, 443 
Buchert T., Dommguez A., Perez-Mercader J., 1999, A&A, 349, 
343 

Burgers J. M., 1948, Adv. Appl. Mech., 1, 171 
Burgers J. M., 1974, The nonlinear diffusion equation : asymp¬ 
totic solutions and statistical problems. D. Reidel Pub. Co. 
cl974, Dordrecht-Holland, Boston, first published in 1973 under 
title: Statistical problems connected with asymptotic solutions of 
the one-dimensional nonlinear diffusion equation 
Cautun M., van de Weygaert R., Jones B. J. T., 2013, MNRAS, 
429, 1286 

Cautun M., van de Weygaert R., Jones B. J. T., Frenk C. S., 2014, 
MNRAS, 441, 2923 

Codis S., Pichon C., Devriendt J., Slyz A., Pogosyan D., Dubois 
Y, Sousbie T., 2012, MNRAS, 427, 3320 
Colberg J. M., Krughoff K. S., Connolly A. J., 2005, MNRAS, 
359, 272 

Coles P, Melott A. L., Shandarin S. R, 1993, MNRAS, 260, 765 
Colless M., Dalton G., Maddox S., Sutherland W., Norberg P, 
et ah, 2001, MNRAS, 328, 1039 

Dietrich J. R, Werner N., Clowe D., Finoguenov A., Kitching T, 
Miller L., Simionescu A., 2012, Nature, 487, 202 
Dolag K., Meneghetti M., Moscardini L., Rasia E., Bonaldi A., 
2006, MNRAS, 370, 656 
Dommguez A., 2000, Phys. Rev. D, 62, 103501 
Dommguez-Tenreiro R., Onorbe J., Martlnez-Serrano F, Serna 
A., 2011, MNRAS, 413, 3022 

Doroshkevich A. G., Ryaben’kii V. S., Shandarin S. R, 1973, As¬ 
trophysics, 9, 144 

Dubois Y. et ah, 2014, MNRAS, 444, 1453 
Dunkley J. et al., 2009, ApJS, 180, 306 

Rorero-Romero J. E., Hoffman Y, Gottldber S., Klypin A., Yepes 
G., 2009, MNRAS, 396, 1815 

Franx M., Illingworth G., de Zeeuw T., 1991, ApJ, 383, 112 

Geller M. J., Huchra J. P, 1989, Science, 246, 897 

Genovese C. R., Perone-Pacifico M., Verdinelli I., Wasserman L., 































22 S. Robles, R. Dommguez-Tenreiro, J. Onorbe and F. J. Martmez-Serrano 


2012, Journal of the American Statistical Association, 107, 788 
Gerhard O. E., 1983, MNRAS, 202, 1159 
Godlowski W., Panko E., Flin R, 2011, Acta Physica Polonica B, 
42, 2313 

Gonzalez R. E., Padilla N. D., 2010, MNRAS, 407, 1449 
Gonzalez-Garcia A. C., Onorbe J., Dommguez-Tenreiro R., 
Gomez-Flechoso M. A., 2009, A&A, 497, 35 
Gonzalez-Garcia A. C., van Albada T. S., 2005, MNRAS, 361, 
1030 

Gott, III J. R., Juric M., Schlegel D., Hoyle R, Vogeley M., 
Tegmark M., Bahcall N., Brinkmann J., 2005, ApJ, 624, 463 
Gurbatov S. N., Malakhov A. N., Saichev A. I., 1992, Nonlinear 
Random Waves and Turbulence in Nondispersive Media (Non¬ 
linear Science: Theory & Applications). John Wiley and Sons 
Ltd 

Gurbatov S. N., Saichev A. I., 1984, Radiofizika, 27, 456 
Gurbatov S. N., Saichev A. I., Shandarin S. R, 1989, MNRAS, 
236, 385 

Gurbatov S. N., Saichev A. I., Shandarin S. R, 2012, Physics Us- 
pekhi, 55, 223 

Hahn O., Carollo C. M., Porciani C., Dekel A., 2007a, MNRAS, 
381,41 

Hahn O., Porciani C., Carollo C. M., Dekel A., 2007b, MNRAS, 
375, 489 

Hahn O., Porciani C., Dekel A., Carollo C. M., 2009, MNRAS, 
398, 1742 

Hidding J., Shandarin S. R, van de Weygaert R., 2014, MNRAS, 
437, 3442 

Hoffman Y., Metuki O., Yepes G., Gottlbber S., Rorero-Romero 
J. E., Libeskind N. I., Knebe A., 2012, MNRAS, 425, 2049 
Huchra J. et ah, 2005, in Astronomical Society of the Pacific Con¬ 
ference Series, Vol. 329, Nearby Large-Scale Structures and the 
Zone of Avoidance, Fairall A. R, Woudt P. A., eds., p. 135 
IckeV., 1973, A&A, 27, 1 
Jenkins A. et ah, 1998, ApJ, 499, 20 

Jones B. J. T., van de Weygaert R., Aragon-Calvo M. A., 2010, 
MNRAS, 408, 897 

Jones D. H., Saunders W., Colless M., Read M. A., Parker Q. A., 
et ak, 2004, MNRAS, 355, 747 

Kaastra J., Finoguenov A., Nicastro R, Branchini E., et ak, 2013, 
ArXiv e-prints 

Kartaltepe J. S. et ak, 2007, ApJS, 172, 320 
Kauffmann G. et ak, 2003, MNRAS, 341, 54 
Kimm T., Devriendt J., Slyz A., Pichon C., Kassin S. A., Dubois 
Y, 2011, ArXiv e-prints 

Kirshner R. P, Oemler, Jr. A., Schechter P. L., Shectman S. A., 
1981, ApJ, 248, L57 

Kirshner R. P, Oemler, Jr. A., Schechter P. L., Shectman S. A., 
1987, ApJ, 314, 493 

Kitzbichler M. G., White S. D. M., 2008, MNRAS, 391, 1489 
Kofman L., Pogosyan D., Shandarin S., 1990, MNRAS, 242, 200 
Libeskind N. I., Hoffman Y, Forero-Romero J., Gottlbber S., 
Knebe A., Steinmetz M., Klypin A., 2013, MNRAS, 428, 2489 
Libeskind N. I., Hoffman Y, Knebe A., Steinmetz M., Gottlbber 
S., Metuki O., Yepes G., 2012, MNRAS, 421, L137 
Lin C. C., Mestel L., Shu F. H., 1965, ApJ, 142, 1431 
Lopez-Sanjuan C. et ak, 2015, A&A, 576, A53 
Lotz J. M., Jonsson P, Cox T. J., Primack J. R., 2010, MNRAS, 
404, 575 

Martmez-Serrano F. J., Serna A., Dommguez-Tenreiro R., Molla 
M., 2008, MNRAS, 388, 39 

Massey R., Rhodes J., Ellis R., Scoville N., Leauthaud A., et ak. 


2007, Nature, 445, 286 

Melon A. L., Buchert T., Weib A. G., 1995, A&A, 294, 345 
Melott A. L., Pellman T. R, Shandarin S. R, 1994, MNRAS, 269, 
626 

Melott A. L., Shandarin S. R, Weinberg D. H., 1994, ApJ, 428, 28 
Metuki O., Libeskind N. I., Hoffman Y, Crain R. A., Theuns T., 
2015, MNRAS, 446, 1458 

Novikov D., Colombi S., Dore O., 2006, MNRAS, 366, 1201 
Onorbe J., Martmez-Serrano F. J., Dommguez-Tenreiro R., Knebe 
A., Serna A., 2011, ApJ, 732, L32 
Patton D. R., Atfield J. E., 2008, ApJ, 685, 235 
Patton D. R., Carlberg R. G., Marzke R. O., Pritchet C. J., da Costa 
L. N., Pellegrini P. S., 2000, ApJ, 536, 153 
Paz D. J., Stasyszyn R, Padilla N. D., 2008, MNRAS, 389, 1127 
Peebles P. J. E., 1980, The large-scale structure of the universe 
Pichon C., Pogosyan D., Kimm T., Slyz A., Devriendt J., Dubois 
Y, 2011, MNRAS, 418, 2493 

Platen E., van de Weygaert R., Jones B. J. T., 2007, MNRAS, 380, 
551 

Pogosyan D., Bond J. R., Kofman L., Wadsley J., 1998, in Wide 
Field Surveys in Cosmology, Colombi S., Mellier Y, Raban B., 
eds.. Editions Frontieres, Gif-sur-Yvette, France, p. 61 
Porciani C., Dekel A., Hoffman Y, 2002a, MNRAS, 332, 325 
Porciani C., Dekel A., Hoffman Y, 2002b, MNRAS, 332, 339 
Powell L. C., Bournaud R, Chapon D., Devriendt J., Gaibler V., 
Khochfar S., Slyz A., TeyssierR., 2013, in lAU Symposium, Vol. 
295, lAU Symposium, Thomas D., Pasquali A., Ferreras I., eds., 
pp. 13-16 

Powell L. C., Slyz A., Devriendt J., 2011, MNRAS, 414, 3671 
Robaina A. R., Bell E. R, van der Wei A., Somerville R. S., Skel¬ 
ton R. E., McIntosh D. H., Meisenheimer K., Wolf C., 2010, ApJ, 
719, 844 

Rykoff E. S. et ak, 2014, ApJ, 785, 104 
Sahni V., Coles P, 1995, Physics Reports, 262, 1 
Serna A., Dommguez-Tenreiro R., Saiz A., 2003, ApJ, 597, 878 
Shandarin S. R, Zeldovich Y. B., 1989, Reviews of Modern 
Physics, 61, 185 

Shen J., Abel T., Mo H. J., Sheth R. K., 2006, ApJ, 645, 783 
Sheth J. V., 2004, MNRAS, 354, 332 
Sheth R. K., van de Weygaert R., 2004, MNRAS, 350, 517 
Sousbie T., 2011, MNRAS, 414, 350 
Sousbie T., Colombi S., Pichon C., 2009, MNRAS, 393, 457 
Sousbie T., Pichon C., Colombi S., Novikov D., Pogosyan D., 
2008a, MNRAS, 383, 1655 

Sousbie T., Pichon C., Courtois H., Colombi S., Novikov D., 
2008b, ApJ, 672, LI 

Sousbie T., Pichon C., Kawahara H., 2011, MNRAS, 414, 384 
Springel V., White S. D. M., Hemquist L., 2004, in lAU Sympo¬ 
sium, Vol. 220, Dark Matter in Galaxies, Ryder S., Pisano D., 
Walker M., Freeman K., eds., p. 421 
Springel V. et ak, 2005, Nature, 435, 629 
Stoica R. S., Martinez V. J., Mateu J., Saar E., 2005, A&A, 434, 
423 

Stoica R. S., Martinez V. J., Saar E., 2007, Journal of the Royal 
Statistical Society: Series C (Applied Statistics) 56 (4), 459-477, 
56, 1 

Stoica R. S., Martinez V. J., Saar E., 2010, A&A, 510, A38 
Tasca L. A. M., Le Fevre O., Lopez-Sanjuan C., Wang P.-W., 
etak, 2014, A&A, 565, AlO 

Tegmark M., Blanton M. R., Strauss M. A., Hoyle R, SDSS Col¬ 
laboration, 2004, ApJ, 606, 702 
Tejos N., 2014, ArXiv e-prints 



Lagrangian Volume Deformations 23 


Tempel E., Stoica R. S., Martmez V. J., Liivamagi L. J., Castellan 

G. , Saar E., 2014, MNRAS, 438, 3465 

Tillson H., Devriendt J., Slyz A., Miller L., Pichon C., 2015, MN¬ 
RAS, 449, 4363 

van de Weygaert R., Bond J. R., 2008, in Lecture Notes in Physics, 
Berlin Springer Verlag, Vol. 740, A Pan-Chromatic View of 
Clusters of Galaxies and the Large-Scale Structure, Plionis M., 
Lopez-Cruz O., Hughes D., eds., p. 335 
van Haarlem M., van de Weygaert R., 1993, ApJ, 418, 544 
Van Waerbeke L., Benjamin J., Erben T., Heymans C., Hilde¬ 
brand! H., et ah, 2013, MNRAS, 433, 3373 
Vergassola M., Dubrulle B., Frisch U., Noullez A., 1994, A&A, 
289, 325 

Weinberg D. H., Hernquist L., Katz N., 1997, ApJ, 477, 8 
Werner N., Finoguenov A., Kaastra J. S., Simionescu A., Dietrich 
J. P, Vink J., Bohringer H., 2008, A&A, 482, L29 
White S. D. M., Silk J., 1979, ApJ, 231, 1 
Wu Y., Batuski D. J., Khalil A., 2009, ApJ, 707, 1160 
Yepes G., Dommguez-Tenreiro R., Couchman H. M. R, 1992, 
ApJ, 401, 40 

Yoshisato A., Matsubara T., Morikawa M., 1998, ApJ, 498, 48 
Yoshisato A., Morikawa M., Gouda N., Mouri H., 2006, ApJ, 637, 
555 

Zel’dovich Y. B., 1970, A&A, 5, 84 

Zhang Y, Yang X., Faltenbacher A., Springel V, Lin W., Wang 

H. , 2009, ApJ, 706, 747 



